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Abstract 


This  work  presents  a  new  approach  for  building  the  feedback  solution 
for  the  minimum  delay  dynamic  message  routing  problem  for  single  destination 
networks . 

The  approach  fully  exploits  the  special  structure  of  the  constraint 
matrices  obtained  in  the  dynamic  state  space  model  suggested  in  previous  works, 
by  transforming  every  linear  program  arising  from  the  necessary  conditions, 
into  a  maximal  weighted  flow  problem. 

Taking  advantage  of  several  properties  concerning  the  networks  cor¬ 
responding  to  the  linear  programming  problems,  all  theorems  regarding  certain 
simplifying  characteristics  of  the  feedback  solution  that  apply  in  the  case 
of  single-destination  networks  with  all  unity  weightings  in  the  cost  functional 
are  reproved  in  a  simplified  and  more  straightforward  manner. 

A  compact  algorithm  for  the  construction  of  the  feedback  solution  is 
presented,  the  algorithm  being  implementable  on  networks  of  reasonable  size. 

A  method  for  obtaining  all  solutions  of  the  linear  programming  prob¬ 
lems  required  by  the  algorithm,  based  on  the  application  of  linear  programming 
techniques  in  networks  is  provided.  The  method  is  implemented  by  a  computer 
program  and  several  examples  are  run  to  test  its  applicability. 

In  addition,  a  deep  geometrical  insight  to  every  step  of  the  algor¬ 
ithm  is  given  by  deriving  the  explicit  set  of  inequalities  defining  the  prob¬ 
lem  constraint  figure  in  the  state-velocity  space. 

The  complexity  of  the  problem  is  also  analyzed,  being  exponential  in 
the  number  of  the  network  nodes,  thus  giving  an  idea  of  the  maximal  network 
size  for  which  a  full  feedback  solution  can  be  obtained  under  the  available 
computational  resources. 


Glossary  of  Notations 


Notation 


Definition 


(i,k) 


P 

(X.X) 


set  of  network  nodes  not  including  the  destination  node 

set  of  network  links 

directed  link  from  node  i  to  node  k 

capacity  of  link  (i,k) 

vector  of  state  variables 

vector  of  velocity  of  state  variables  C-x) 

vector  of  control  variables 

vector  of  costates 

constraint  figure  in  u-space 

constraint  figure  in  y- space 

set  of  states  travelling  on  interior  arcs  on  [tp,tp+1) 
set  of  states  travelling  on  boundary  arcs  on  [t  ,t  j.) 
set  of  states  leaving  the  boundary  at  tp  backwards  in  time 
feedback  control  region  constructed  from  optimal  traject¬ 
ories  on  the  segment  [tp,tp^) 
set  of  operating  points  on  [t  ,tp+p 
a  minimal  cut  of  the  network 
demand  for  flow  in  node  i 


supply  of  flow  to  the  network 
cardinality  of  I 
cardinality  of  Lp 
the  Hamiltonian  at  x 
Convex  Hull 

Collection  of  nodes  k  such  that  (i,k)  e  L 
Collection  of  nodes  l  such  that  (A,i)  c  L 
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Section  1 

INTRODUCTION 

In  [SlJ  a  state  space  model  for  dynamic  routing  in  data  communication 
networks  is  suggested.  The  main  feature  of  this  model  is  that  it  permits  to 
express  the  delay  experienced  by  the  messages  travelling  in  the  network  in 
terms  of  state  and  control  variables  describing  the  problem  instead  of  models 
based  on  queueing  theory.  The  latter  requires  explicit  closed-form  expres¬ 
sions  for  the  average  delays  which  can  be  found  analytically  only  for  very 
special  distributions  and  dependence  relationships.  The  model  also  permits 
to  develop  closed  loop  strategies  for  the  message  routing  problem  and  can 
handle  transients  by  changing  the  routing  policy  in  a  dynamic  fashion. 

We  begin  by  presenting  a  brief  description  of  the  model,  simplified 
to  the  case  of  single  destination  networks  with  zero  inputs.  Consider  first 
the  following  notations: 

N  =  {l,2,...,n}  is  the  set  of  network  nodes  (not  including  the  destination  node), 
d  =  destination  node. 

L  =  { (i,k)/i,k  e  WU  d  and  there  is  a  direct  link  connecting  i  to  k}  is  the 

set  of  network  links, 

E(i)  =  collection  of  nodes  k  such  that  (i,k)  e  L. 

I(i)  *  collection  of  nodes  l  such  that  (i,i)  e  L. 

All  links  of  the  network  are  taken  to  be  unidirectional.  Now,  looking  at  the 
network  from  a  macroscopic  point  of  view,  the  number  of  messages  in  each  node 
can  be  approximated  by  a  continuous  variable  called  "amount  of  traffic".  The 
state  variables  of  the  system  are  defined  as  follows: 
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x^(t)  ■  amount  of  traffic  at  node  i,  at  time  t,  where  ieW. 

The  control  variables  are  defined  as: 

uik(t)  »  amount  of  traffic  on  link  (i,k)  at  time  t  where  (i,k)  el. 
The  dynamics  of  the  system  are  given  by  the  equations  describing  the  rate  of 


change  of  the  contents  of  each  node,  namely: 

x  (t)  =  -  l  u  (t)  +  l  u  (t) 
1  keE(i)  1K  lei  (i) 


y  ieW 


The  constraints  are: 


x.(t)  >  0  , 


(1.1) 


(1.2a) 


where 


U:  0  <  u.  (t)  <  c._  , 

-  m  —  lk 


(1.2b) 


c^  *  capacity  of  link  (i,k)  in  units  of  traffic/unit  time,  where 
(i.k)  e  L. 

The  cost  functional  is  taken  to  be  the  total  delay  experienced  by  the 
messages  travelling  in  the  network,  starting  at  a  given  time  t  and  ending  at 
a  time  t^  when  the  network  is  emptied,  i  e  x^(t-)  *  0,  Vic  W.  The  ex¬ 
pression  for  the  above  quantity  is 


a. 

f  [  l  xi (t) ]dt  . 

i  ieN  1 


(1.3) 


From  now  on,  the  various  variables  of  the  model  will  be  represented  in  vector 
form  as  follows:  Denoting  by  x  and  u  the  respective  concatenations  of 
state  variables  and  controls,  to  every  differential  constraint  in  (1.1)  cor¬ 
responds  a  vector  such  that: 
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where  the  vectors  are  the  rows  of  the  incidence  matrix  B  of  the  network. 

Now  we  express  the  linear  optimal  control  problem  with  linear  state 
and  control  variable  inequality  constraints  representing  the  data  communication 
network  closed  loop  dynamic  routing  problem  as  stated  in  [SI]: 

Find  the  set  of  controls  u  as  a  function  of  time  and  state, 
u(t)  *  u(t,x)  ,  t  e  '[tQ,tf]  , 

that  brings  a  given  initial  condition  x(tQ)  *  to  the  final  condition 
x(t^)  =  0  and  minimizes  the  cost  functional  (1.3]  subject  to  the  dynamics  (1.1) 
and  to  the  state  and  control  variable  inequality  constraints  (1.2). 

It  is  well  know  that  in  most  optimal  control  problems  it  is  quite 
difficult  to  obtain  feedback  solutions.  In  our  case  however,  owing  to  the 
linearity  of  the  problem,  we’ can  cover  the  entire  state  space  with  optimal 
controls  by  solving  just  a  comprehensive  set  of  linear  programs.  In  [SI] 
an  approach  is  suggested,  by  way  of  a  simple  example  for  constructing  the 
feedback  solution.  In  [Ml],  [M2]  and  [M3]  this  approach  is  elaborated  upon 
by  developing  the  so  called  Constructive  Dynamic  Programming  Algorithm  for 
the  construction  of  the  feedback  solution.  We  proceed  now  with  a  brief  pre¬ 
sentation  of  the  results  obtained  in  the  above  works  and  the  conceptual 
lines  of  the  Constructive  Dynamic  Programming  Algorithm  in  order  to  provide 
the  reader  with  basic  notions  for  the  understanding  of  the  present  report. 

We  begin  by  presenting  the  necessary  and  sufficient  conditions  of  optimality: 

Theorem  1.1:  (See  [Ml,  pp.  53-63]  or  [M2,  pp.  13-20]). 

Let  the  scalar  functional  h  be  defined  as  follows: 

h(u(t),X(t))  £  XT (t)x(t)  *  XT(t)Bu(t)  . 
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A  necessary  and  sufficient  condition  for  the  control  law  u*(*)  £  U  to  be 
optimal  for  problem  (1.1)  -  (1.3)  is  that  it  minimizes  h  pointwise  in  time 
namely: 

AT(t)Bu*(t)  <,  AT(t)Bu(t)  (1.4) 

V  u(t)  £  U  ,  V  t  c  [tQ,tf]  . 

The  costate  A(t)  is  possibly  a  discontinuous  function  which  satisfies  the 
following  differential  equations: 

-dAi(t)  =dt  +  dui(t),  VieW  tE[tQ,tf]  ,  (1.5) 

where  componentwise  du^(t)  satisfies  the  complementary  slackness  conditions 

(1.6) 

The  terminal  boundary  condition  for  the  costate  differential  equation  is 

A(t^)  =  0  free  ,  (1.7) 

and  the  transversality  condition  is 

XT(tf)x(tf)  =  0  .  (1.8) 

Finally,  the  function  h  is  everywhere  continuous,  i.e.: 

h(u(t"),A(t"))  =  h(u(t+) ,  A(t+) ) ,  Vte[tQ,tf]  . 

□  Theorem  1 . 1 

From  inequality  (1.4)  of  the  necessary  and  sufficient  conditions  we 
see  that  the  optimal  control  function  u*(*)  is  given  at  every  time 
t  e  by  the  solution  of  the  linear  program 

u*(x)  *  ARGMIN[AT(t)BU'(t)  ] 
u(x)etl 


xi(t)dui(t)  =  0 


du.  (t)  <  0 


V  t  £  [t0,tf] 

V  i  e  N 


(1.9) 
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From  (1.9),  and  owing  to  the  structure  of  the  incidence  matrix  B,  it  follows 
that  there  always  exists  an  optimal  solution  for  which  the  controls  are  piece- 
wise  constant  in  time.  Moreover,  since  the  equations  governing  the  dynamics 
of  the  system  are  linear,  the  corresponding  state  trajectories  have  piecewise 
constant  slopes.  From  the  nature  of  the  controls  follows  that  every  optimal 
trajectory  may  be  characterized  by  a  finite  number  of  parameters.  These  para¬ 
meters  are: 

U(x)  =  (u  , u. , . . . , u «  _ } 

-o  -1  r-1 

and 

T(x)  -  (to,t1,...,tf>  , 

where  U(x)  is  the  sequence  of  optimal  controls,  T(x)  is  its  associated 
control  switch  time  sequence,  and  the  element  is  the  optimal  control  on 

"[yVii-  . f-‘i- 

Moreover,  every  segment  of  an  optimal  trajectory,  say  the  segment 
[t  ,tp+j]  is  characterized  by  the  following  parameters: 

BP  =  ui/xi(T)  s  0  VTeIWi)} 

and 

!P  4  «vxi(T)  * 0  VTttyyi)>* 

that  is.  Bp  and  1^  are  the  sets  of  states  travelling  on  boundary  arcs  and 
interior  arcs  respectively  on  [tp,t  j).  Now,  the  main  fact  following  from 
the  nature  of  the  optimal  trajectories  is  that  the  state  space  can  be  divided 
into  regions,  each  one  being  a  convex  polyhedral  cone,  when  to  every  point  of 
a  specific  region  corresponds  an  identical  set  of  optimal  controls.  This  is 
exactly  the  reason  that  makes  the  construction  of  a  feedback  solution  possible 
The  above  regions  are  referred  to  as  "feedback  control  regions"  and  are  de¬ 
noted  by  R. 

The  problem  is  then  to  construct  these  feedback  control  regions,  there 
by  specifying  the  optimal  control  for  every  point  of  the  state  space.  Now, 
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a 


considering  the  special  geometric  characterization  of  the  feedback  space,  and 
thinking  in  the  spirit  of  dynamic  programming,  let  us  look  at  the  optimal 
trajectories  backwards  in  time,  beginning  at  tf.  We  will  then  see  a  sequence 
of  states  leaving  the  boundary  (perhaps  two  or  more  at  a  time)  and  varying 
with  constant  slopes.  Now,  in  [M2]  it  is  proven  that  if  any  state  variable, 
say  is  strictly  positive  on  the  last  time  interval  [t^  j.t^]  of  an 

optimal  trajectory,  then  X^(tp  *  0.  Moreover,  by  (1.6)  we  know  that  for 
this  state  we  have  dvu(r)  =  0  Vie  [t^  ^,t^]  50  that  by  (1.5)  we  obtain 

Xj,(T)  *  -1  V  t  e  [tfl,tf] . 

Now  suppose  we  had  decided  to  check  if  there  exist  optimal  traject¬ 
ories  in  the  state  space  with  a  specified  set  of  states  travelling  on 

interior  arcs  in  the  last  interval  and  with  the  remaining  states  travelling 
on  boundary  arcs.  Since  we  know  the  costates  corresponding  to  the  states  in 
If  A  (by  the  above  arguments)  and  using  (1.9),  these  trajectories  may  be 
found  by  solving  the  following  linear  program: 


s.t. 


Find  all: 


u*  =  ARGMIN  l  X.(x)i  a  ARGMIN  l  X.  (r)b.u 

uell  X.el-  ,  1  1  ueU  x.el.  .  1  ~1- 

i  f-1  -  l  f-1 


x.  <  0 

l  — 

Vx.elf.j 

O 

If 

•  >< 

V  x.  e  B-  . 

J 

i  f-1 

T  e 

(1.10) 


The  linear  program  (1.10)  is  called  the  constrained  optimization  prob¬ 
lem  and  its  solutions  are  trajectories  only  provided  that  there  exist  values 
Xj  VXj  that  satisfy  the  necessary  conditions  (1.5)  -  (1.6).  Moreover, 

the  solutions  of  (1.10)  provide  optimal  directions  xT  in 

state  space,  defining  the  convex  polyhedral  cone  corresponding  to  the  feedback 
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control  region  characterized  by  the  set  of  controls  that  solves  (1.10).  Now, 
solving  such  a  linear  program  for  every  possible  combination  of  state  vari¬ 
ables  leaving  the  boundary  (backwards  in  time)  at  t^,  we  obtain  feedback 
control  regions  containing  all  points  from  which  the  origin  can  be  reached 
while  maintaining  sets  1^  ^  away  from  the  boundary  for  all  t  <  t^.  In  fact 
the  linear  program  (1.10)  must  be  solved  parametrically  in  time  until  the  con¬ 
trol  changes.  A  change  in  the  control  defines  a  hyperplane  passing  through 
the  origin  called  "breakwall".  By  solving  (1.10)  until  the  control  does  not 
break  any  more,  we  obtain  a  convex  polyhedral  cone  divided  by  breakwalls  into 
a  finite  number  of  regions.  Each  region  is  characterized  by  a  specific  control 
set  and  is  referred  to  as  "break  feedback  control  region".  The  last  region 
constructed  when  solving  (1.10),  (i.e.  when  there  are  no  more  breaks  in  the 
control)  is  called  "non-break  feedback  control  region". 

In  order  to  continue  with  the  description  of  the  algorithm,  let  us 
introduce  the  following  definition: 

L  =  (x^/x^  g Bp  and  x^  is  designated  to  leave  the  boundary 

backward  in  time  at  t  }  . 

P 

Now,  in  order  to  cover  the  whole  state  space  with  optimal  controls, 
we  must  take  every  feedback  control  region  already  constructed  and  allow  every 
possible  combination  of  states  still  travelling  on  the  boundary  to  leave  it, 
backward  in  time.  In  general,  this  step  is  performed  as  follows:  Denoting 
by  Rp  the  feedback  control  region  constructed  from  the  optimal  trajectories 
on  the  segment  [tp,tpl),  pick  a  set  of  states  of  Bp,  say  Lp  and  parti¬ 
tion  Rp  into  "subregions"  with  respect  to  Lp.  The  concept  "subregion" 
will  be  explained  later.  In  order  to  let  Lp  leave  the  boundary  at  a  certain 

time  t  (called  boundary  junction  time),  we  must  find  the  costate  values 
P 

X,  (t  )  Vx,  e  L  that  allow  the  departure  of  the  states  in  L  from  the 
i  p  i  p  r  p 
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boundary  while  still  being  optimal.  The  above  costate  values  are  referred  to 
as  "leave- the-boundary  costate  values".  In  geometrical  terms  the  above  co¬ 
state  values  are  found  as  follows: 


By  defining 


y  =  -x  =  -Bu 


y  -  {ye  Rn/u  e  U>  , 


we  can  transform  the  linear  program  {1.9)  into  the  following  linear  program 
with  decision  vector  y{t): 


y*(x)  *  ARGMAX  X1(t)v(t) 

y(f)ey 


(1.11) 


so  that  at  t^  we  can  consider  the  hyperplane  given  by  z(t^)  =  X  (t  )y 
(called  the  Hamiltonian)  tangent  to  V  (called  the  y-constraint  figure)  that 
provides  points  of  tangency  (called  operational  points)  that  are  in  fact 

the  optimal  directions  in  the  state  space  defining  the  feedback  control 
region  R^.  Now  in  order  to  find  the  leave -the -boundary  costates  for  L^, 
we  rotate  the  Hamiltonian  around  until  we  touch  a  surface  of  tangency 
of  y  called  l^-positive  face,  having  at  least  one  point  with 

y.  >  0  V  x.  e  L  . 

7\  l  p 

The  new  orientation  of  the  Hamiltonian  gives  the  desired  leave-the-boundary 
costates  and  now  We  solve  a  linear  program  of  the  form: 

Find  all: 


u*  -  ARGMIN  I  X. (t)x.  *  ARGMIN  £  X. (r)b.u 
_ a  i  i  -  _t  i 


x.  el  . 

l  p-1 


x.  el  , 

l  p-1 


x.  <  0  V  x.  e  I  .  »  I  UL 
i  —  l  p-1  p  t 


(1.12) 


x  «  0  V  x.  e  B  ■  B  /L 
i  i  p-1  p  p 


V  t  e  (t  ,t  ) 

p  p-1 
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where  (1.12)  is  again  solved  parametrically  in  time  until  the  control  breaks. 


A  subregion  of  is  the  set  of  all  points  which,  when  taken  as  the 
point  of  departure  of  L^,  result  in  a  common  set  of  optimal  controls. 
Clearly,  such  a  partitioning  of  may  exist  since  to  two  distinct  points  in 
Rp  correspond  different  leave-the-boundary  costates.  In  geometrical  terms, 
we  could  see  the  Hamiltonian  rotating  according  to  the  change  in  the  costates 


when  the  states  are  travelling  in  R^.  When  allowing  to  leave  the 
boundary  from  distinct  points  of  R^,  the  Hamiltonian  could  meet  different 
Lp-positive  faces,  therefore  giving  different  leave-the-boundary  costate  values. 


The  steps  described  for  the  construction  of  feedback  control  regions 
form  the  so-called  Constructive  Dynamic  Programming  Algorithm  as  stated  in 
[Ml],  [M2]  and  [M3].  Several  features  make  the  algorithm  to  be  a  conceptual 
one  rather  than  an  implementable  one.  Among  them,  breaks  in  the  optimal  con¬ 
trols,  the  existence  of  the  subregions  just  pointed  out,  rion-uniqueness  of  the 
leave-the-boundary  costates,  and  non-global  optimality  of  certain  sequences, 
in  addition  to  computational  complexities  associated  with  the  algorithm,  for 
instance  the  problem  of  finding  all  solutions  to  linear  programs.  However, 
it  turns  out  that  when  dealing  with  single  destination  networks,  and  when  no 
priorities  are  assigned  to  the  nodes  (that  is  the  cost  functional  has  unity 
weightings),  many  simplifications  are  attained,  leading  to  a  compact  and  imple¬ 
mentable  algorithm  at  least  for  moderate  size  networks. 


In  the  present  paper,  we  deal  with  the  construction  of  the  feedback 
solution  to  the  minimum  delay  dynamic  message  routing  problem  for  single 
destination  networks,  with  zero  inputs  and  when  the  cost  functional  has  unity 
weightings.  The  approach  used  to  tackle  the  problem  is  based  on  the  trans- 
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formation  of  every  linear  program  into  a  maximal  weighted  flow  problem.  The 
transformation  permits  not  only  to  develop  a  compact  algorithm  for  obtaining 
the  feedback  solution,  but  also  to  prove  in  a  straightforward  manner  all  the 
simplifying  features  characterizing  the  feedback  space,  to  develop  a  suitable 
method  for  finding  all  the  solutions  to  the  linear  programming  problems  and 
to  obtain  explicitly  the  constraint  figure  in  the  state  velocity  space  for 
the  geometrical  understanding  of  the  algorithm.  The  approach  also  permits 
to  gain  insight  into  the  complexity  of  the  problem,  thus  providing  the 
essential  information  needed  to  estimate  the  maximal  size  of  the  networks  for 
which  a  complete  feedback  solution  can  be  obtained  under  the  available  compu¬ 
tational  resources. 

One  of  the  most  interesting  features  of  the  algorithm  presented  here 
is  that  although  the  efficiency  of  the  method  for  obtaining  all  optimal  solu¬ 
tions  to  the  linear  programs  is  reduced  by  the  high  degeneracy  of  the  problems 
at  hand,  this  is  compensated  by  the  fact  that  the  higher  the  degree  of 

« 

degeneracy,  the  lower  is  the  number  of  linear  programs  that  have  to  be  solved. 

Moreover,  the  algorithm  requires  to  check  all  possible  sequences  of  states 
leaving  the  boundary  backward  in  time  or,  in  other  words,  all  the  possible 
trajectories  in  the  state  space  to  insure  the  complete  covering  of  the  state 
space.  Again  the  higher  the  degree  of  degeneracy,  the  lower  the  number  of 
such  sequences  that  have  to  be  actually  checked,  thus  obtaining  a  further 
reduction  in  the  complexity  of  the  algorithm 

The  organization  of  the  paper  is  as  follows: 

In  Section  2,  a  general  description  of  the  algorithm  is  presented  and 
the  transformation  of  the  linear  programs  into  maximal  weighted  flow  problems 
is  introduced.  Section  3  is  devoted  to  the  theoretical  results  obtained, 
namely,  the  structure  of  the  y-constraint  figure  and  the  proofs  regarding  all 

J 
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simplifying  features  of  the  feedback  solution  that  apply  in  the  case  of  single 
destination  networks. 

In  Section  4  the  algorithm  is  presented  in  a  compact  form,  the  method 
for  obtaining  all  solutions  to  the  linear  programs  is  described  and  an  example 
of  the  construction  of  the  feedback  solution  is  given. 

Finally,  in  Section  S  a  brief  discussion  about  the  contributions  of 
this  work  is  carried  out  and  topics  for  further  research  in  the  area  are 
suggested. 
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Section  2 


GENERAL  DESCRIPTION  OF  THE  ALGORITHM 

We  begin  with  a  brief  preview  of  this  section.  In  part  A  we  present 
a  simplified  mathematical  description  of  the  Constructive  Dynamic  Programming 
Algorithm  for  building  the  feedback  solution  for  the  minimum  delay  dynamic 
message  routing  problem  for  single  destination  networks  with  all  unity  weight 
ings  in  the  cost  functional.  The  algorithm  consists  of  two  steps.  In  the 
first  step,  we  construct  feedback  control  regions  resulting  from  the  states 
leaving  the  boundary  at  the  final  time  (backward  in  time) .  The  method  to  con 
struct  these  regions  is  depicted  in  parts  B  and  C  of  this  section.  In  the 
second  step,  the  rest  of  the  state  space  is  filled  with  optimal  controls  by 
constructing  the  feedback  control  regions  resulting  from  states  leaving  the 
boundary,  starting  from  already  constructed  regions.  Part  D  of  this  section 
deals  with  the  construction  of  the  above  regions. 

A.  Mathematical  Statement  of  the  Algorithm 

In  [M2],  the  following  properties  associated  with  the  Constructive 
Dynamic  Programming  Algorithm  are  mentioned: 

(a)  Non-global  optimality  of  certain  sequences  of  state  variables  leaving 
the  boundary  backwards  in  time. 

(b)  Non-uniqueness  of  leave- the-boundary  costates. 

(c)  Subregions. 

(d)  Return  of  state  variables  to  boundary  backwards  in  time. 

(e)  Breaks  in  the  optimal  control  between  boundary  junction  times. 
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In  [M3],  a  discussion  of  the  above  properties  is  presented  in  geo¬ 
metrical  terms.  These  problems  complicate  the  formulation  of  a  computational 
scheme  to  implement  the  algorithm.  However,  as  we  will  prove  in  Section  3  of 
this  paper,  it  turns  out  that  these  problems  do  not  apply  in  the  case  of 
single  destination  networks  with  all  unity  weightings  in  the  cost  functional, 
thus  permitting  the  following  simplified  statement  of  the  algorithm: 

Step  1:  Solve  the  following  series  of  linear  programming  problems: 

Find  all 

s.t. 

for  all 
when  W 


u*  =  ARGMIN  l  X.  =  ARGMIN  £  b.u 
ueli  x.eLi.  1  u  e  U  x.eL* 


x.  <  0  Vx.eL- 
i  —  if 


x.  =  0  V  x.eB 
3  3  f-1 


lf  C  {xVieN} 

denotes  the  set  of  nodes  not  including  the  destination. 


Step  2:  For  all  feedback  control  regions  R^: 

(a)  Calculate  the  leave-the-boundary  costate  values  X^(tp)  for  all 
x.^  e  Bp,  tp  being  an  arbitrary  boundary  junction  time. 

(b)  Solve  the  following  series  of  linear  programming  problems: 

Find  all: 


u*  *  ARGMIN  I  X.(t  ) x.  =  ARGMIN  T  X.  (t  )b  u 

u  e  U  x.cl  ,  1  P  1  »»U  X.  El  ,  1  P  -1' 

1  P-1  -  1  P-1 

x.  <  0  V  X.  e  I  , 

l  —  l  p-1 

x.  *  0  V  x.  e  B 

i  i  P-1 


(2.2) 


for  all 


IS 


As  we  said  before,  in  the  first  step,  feedback  control  regions  result¬ 
ing  from  the  states  leaving  the  boundary  at  t^,  are  constructed.  Note  that 
in  problem  (2.1),  the  cost  function  has  unity  weightings.  The  reason  is  that 
by  Corollary  2  in  [M2],  A^(t^)  *  0  V e  L^,  and  by  the  necessary  and  suffici¬ 
ent  conditions  for  optimality  (see  [M2,  pp.  13-17]),  A^(x)  *  -1  Vx^ e  Lf. 

Moreover,  since  there  are  no  breaks  between  boundary  junction  times,  the  co¬ 
states  corresponding  to  states  in  are  equal  for  all  tc  (-*,tf)  so  that 
at  every  time  in  this  interval  the  cost  function  is  the  same  after  normaliza¬ 
tion.  From  this  fact  follows  that  at  every  time  in  (-*,t^)  we  obtain  the 
same  set  of  optimal  solutions  to  problem  (2.1),  so  that  the  time  plays  no  role 
here.  Notice  also,  that  since  each  set  of  state  variables  leaving  the  boundary 
corresponds  to  a  globally  optimal  solution,  every  solution  to  (2.1)  gives  an 
optimal  trajectory. 

In  the  second  step,  feedback  control  regions  resulting  from  states 
leaving  the  boundary  from  previously  constructed  regions  are  built.  Since 
there  is  only  one  subregion  per  region,  it  follows  that  from  every  point  of 
a  given  feedback  control  region,  states  leave  the  boundary  with  the  same  con¬ 
trol  set  and  therefore  boundary  junction  times  may  be  chosen  arbitrarily.  The 
above  reasoning  together  with  the  fact  that  the  leave- the-boundary  costates 
are  unique  for  a  given  boundary  junction  time,  and  with  the  fact  that  no  state 
variable  is  ever  required  by  optimality  to  increase  forward  in  time,  justify 
the  simplifications  obtained  in  the  statement  of  the  algorithm. 

The  first  step  is  carried  out  as  follows:  First,  the  linear  program 
corresponding  to  the  leave-the-boundary  of  all  states  is  solved.  Second,  the 
linear  programs  corresponding  to  (n-1)  states  leaving  the  boundary  are  solved 
(when  n  is  the  number  of  nodes  in  the  network  not  including  the  destination 
node),  and  so  on.  As  it  will  be  shown  later,  this  way  of  solving  the  first 
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step  enables  in  general  a  considerable  reduction  in  the  number  of  linear  pro¬ 
grams  that  have  to  be  solved.  Moreover,  it  will  also  be  shown  later  that  in 
fact  the  second  step  does  not  require  solving  any  linear  programs,  thus  greatly 
reducing  the  complexity  of  the  algorithm. 

B.  All  States  Leave  the  Boundary 


According  to  (2.1),  the  linear  programming  problem  we  have  to  solve  for 
this  case  is: 


Find  all: 


u*  =  ARGMIN  l  X.  =  ARGMIN  £  b.u  , 
ueU  i  eN  1  ueti  i  eN  1 


s.t. 


x.  <  0  V  i  e  W  , 

i  ~ 


U  = 


0  <  u..  <  c., 

-  jk  -  jk 


l  V  (j  ,k)  e  L 


(2.3) 


(2.3a) 


(2.3b) 


Adding  slack  variables  y.  ^  0  to  the  differential  constraints  (2.3a),  prob¬ 
lem  (2.3)  is  transformed  into: 


Find  all: 


u*  *  ARGMIN  l  b.u  *  ARGMAX  £  y.  , 
ueU  icN  ~x~  ueU  ieW 


(2.4) 


s.t. 


^  =  0 


V  i  eN  , 


(2.4a) 


y,  >  0 


0  1  Ujk  i  Cjk 


IV(j.k)  eL 


(2.4b) 


Notice  that  in  teras  of  the  variable  "y",  the  minimization  problem  (2.3)  is 


transformed  into  the  maximization  problem  (2.4). 


* 


Suppose  we  are  interested  in  finding  only  one  solution  to  problem  (2.4). 
In  order  to  do  this,  we  transform  the  linear  program  (2.4)  into  the  following 
Maximal  Flow  Problem: 

Maximize  F  =  £  Y-  >  (2.S) 

ieW  1 


s.t. 


“  l  uii  +  l  ,uki +  yi  =  0 

jeE(i)  kel(i)  1  1 


V  i  e  N 


y.  >  0 

'  l  — 


(2.5a) 


ueU 


<  u..  <  c.t 
-  Jk  -  Jk 


V  (j  ,k)  e  L 


(2.5b) 


The  network-  corresponding  to  problem  (2.5)  is  formed  by  adding  to  the 
original  network  a  new  node  called  the  "source",  and  n  links  with  no  capacity 
constraints  connecting  the  source  with  each  node  of  the  network  (except  for  the 
destination  node).  The  slack  variable  "y/'  represents  the  flow  on  the  link 
connecting  the  source  with  node  i.  F  represents  the  total  flow  into  (and  out 
of)  the  network  (see  Fig.  2.1). 

s  *  source  node 
d  »  destination  node 


The  links  connecting 
the  different  nodes 
are  omitted  here  to 
avoid  causing  havoc 
in  the  figure. 

Figure  2.1  -  Network  representing  Problem  (2.5) 
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The  transformation  of  the-  linear  program  into  a  Maximal  Flow  Problem 
is  the  fundamental  step  required  for  the  proofs  of  the  theorems  regarding  the 
simplifying  properties  of  the  feedback  solution.  In  Appendix  A,  we  give  basic 
concepts  of  graph  theory  that  we  use  in  the  development  of  the  algorithm  and 
in  the  proofs  of  the  theorems. 

Now  we  continue  with  the  solutions  of  problem  (2.5).  Our  aim  is  to 
find  one  solution  to  problem  (2.5)  by  means  of  the  Maximal  Flow  Algorithm  and 
then,  without  changing  the  flow  achieved  in  the  first  optimal  solution,  to  per¬ 
form  all  possible  pivots  on  the  network  in  order  to  find  all  the  remaining 
optimal  solutions  of  problem  (2.5).  Notice  that  since  we  are  seeking  extremal 
optimal  solutions,  we  require  that  the  first  one  will  be  an  "extended  basic 
optimal  solution",  (see  Appendix  A).  In  this  way  every  solution  obtained  by 
pivoting  will  also  be  an  extended  basic  optimal  solution,  thus  representing  an 
extremal  point  in  the  optimization  space.  It  turns  out  that  finding  a  first 
"extended  basic  optimal  solution"  to  problem  (2.5)  is  trivial.  From  the  fact 
that  the  links  corresponding  to  the  "y"  variables  do  not  have  capacity  con¬ 
straints  follows  that  there  exists  only  one  minimal  cut  in  the  network  (denoted 
by  (X,X)),  and  it  is  composed  by  all  the  links  connecting  the  nodes  with  the 
destination;  thus: 

X  =  (s,l,2, . . . ,n)  , 

X  =  (d)  . 

Then  we  choose  as  the  first  optimal  solution 

c.d  iel(d) 

0  otherwise 

uid’Cid  =  ltd) 

-  uki  ■  0  j  e  E(i), 

VieN 


the  following  flows: 


k  e  I (i) 


(2.6) 
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and 


l 

iel(d) 


Cid 


(2.6a) 


Clearly,  (2.6)  is  an  "extended  basic  solution"  if  the  variables  y^ 
are  declared  basic  and  the  remaining  variables  are  declared  non-basic.  Notice 
that  if  there  is  no  direct  connection  between  a  node,  say  node  i,  and  the 
destination  node,  then  the  variable  y.  is  basic  with  value  zero.  Now,  since 

i 

networks  have  generally  many  nodes  not  connected"  directly  with  the  destination 
we  are  dealing  here  with  linear  programming  problems  of  a  high  degree  of 
degeneracy.  The  degeneracy  problem  lowers  the  efficiency  of  linear  program¬ 
ming  methods,  however  in  our  case,  the  higher  the  degree  of  degeneracy,  the 
lower  is  the  number  of  linear  programs  we  have  to  solve.  This  subject  will 
be  discussed  in  Part  C  of  this  section. 


Recall  that  after  finding  the  first  "extended  optimal  basic  solution" 
we  are  interested  in  maintaining  the  achieved  flow  to  assure  that  every  solu¬ 
tion  obtained  by  pivoting  operations  on  the  network  will  be  optimal.  But 
maintaining  maximal  flow  in  the  network  implies  that  at  every  solution  the 
flow  on  the  minimal  cuts  will  be  constant  and  maximal;  thus,  before  pivot¬ 
ing  we  can  reduce  the  number  of  variables  of  the  problem  by  transforming  the 
network  to  the  one  depicted  in  Figure  2.2: 


Figure  2.2  -  Network  configuration  before  the  pivoting  operation. 
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In  Figure  2.2,  accounts  for  the  "demand”  of  flow  in  node  i,  and 

bg  is  the  "supply"  of  flow  to  the  network  where: 


b.  = 
1 


c .  , 
id 


10 


i  e  I  (d) 

Vi  e  N 
otherwise 


b  = 
s 


I 

iel(d) 


c 


id 


(2.7) 


Clearly,  the  reduction  in  the  number  of  links  facilitates  the  work  of  finding 
all  optimal  solutions.  After  the  reduction,  the  remaining  optimal  solutions 
are  found  by  means  of  the  algorithm  we  present  in  Section  4. 


C.  A  Subset  of  States  Leaves  the  Boundary 


In  this  case,  we  let  a  subset  Lf  C. (xVi  e  W},  such  that  ^  (x./i  e  W} 

to  leave  the  boundary  at  tf,  while  all  other  states  Bf  ^  remain  on  the 
boundary.  By  (2.1),  the  linear  program  to  be  solved  is: 

Find  all: 


s.t. 


U*  =  ARGMIN  l  X.  *  ARGMIN  [  b.u 

ueU  x.eL„x  uelJ  x.eL  “1_ 
if  -  if 

x.  <  0  V  x.  z  l.  , 
l  -  if 

■  0  v  £  6f-i  • 


(2.8) 


(2.8a) 


Transforming  (2.8)  into  a  Maximal  Flow  Problem,  we  have: 

Maximize  F  *  £  y.  (2.9) 

x.eL.  1 
i  f 


s.t. 


'  2  “ii  +  I  uki  +  yi  3  0  V  xi  e  Lf 

jsE(i)  1J  kel(i)  K1  1  1  r 


I  uij  -  I  ■  0  V  x 


jeE(i)  1J  kel(i) 


e  B 


f-1 


y.  >  0  Vx.  e  L 

i  —  if 


u  e  U  - 


0  <  u.k  <  c.. 

—  3  —  Jk 


[  V  (j,k)  eL 


(2.9a) 


The  network  corresponding  to  problem  (2.9)  is  formed  by  adding  a  new 
node  called  the  "source"  to  the  original  network,  as  in  Part  B  of  this  section. 
The  difference  here  is  that  we  add  links  without  capacity  constraints  connect¬ 
ing  the  source  only  with  those  nodes  that  correspond  to  state  variables  of  L^. 
Nodes  corresponding  to  state  variables  in  Bf  1  will  be  not  connected  with  the 
source.  For  example.  Figure  2.3  shows  a  network  with  five  nodes  corresponding 
to  the  case  of  two  states  leaving  the  boundary: 


Figure  2.3  -  Example  of  network  when  only  some  states  leave  the  boundary. 

Our  aim  is,  as  in  Part  A,  to  find  a  first  "extended  basic  optimal  solu 


tion"  to  problem  (2.9)  and  then  by  pivoting  to  find  the  remaining  solutions. 
But  here  we  cannot  obtain  the  first  optimal  solution  by  inspection  since  it  is 


not  trivial  to  localize  a  minimal  cut  in  the  network.  Moreover,  we  can  obtain 
an  optimal  solution  by  applying  the  Maximal  Flow  Algorithm  to  the  network,  but 
this  solution  may  not  necessarily  be  an  extended  basic  one.  In  order  to  over¬ 
come  this  difficulty  we  present  now  two  theorems  and  a  corollary  that  will 
assist  us  in  developing  a  method  to  solve  the  problem  when  only  some  states 
leave  the  boundary.  In  addition  these  results  will  show  that  in  general  there 
are  many  combinations  of  state  variables  leaving  the  boundary  for  which  it  is 
not  necessary  to  solve  a  linear  programming  problem. 


Theorem  2 . 1 

Consider  the  following  two  linear  programming  (or  Maximal  Flow)  problems 
Problem  I 

Maximize  F  »  [  y. 

1  x.eL'  1 
l  f 


s.t. 


I  u.  .  +  l  u,  .  +  y. 

jeE(i)  1J  kel (i)  *a  1 

1  u.  .  -  £  u,.  =  0 

jcE(i)  13  kel(i) 


0  Vxt  e  L£ 


V  X.  e  B' 
i  f-1 


Problem  II 


y.  >  0  V  x.  e  L' 

7 1  —  if 

u  e  U  . 

Maximize  ^  2  s  [ 

x. eL,.  1 
l  f 


I  u.  .  +  F  u,  .  +  y. 
jeE(i)  13  kel(i)^1  1 

I  ^  *  l  =  0 

jeE(i)  13  kel  (i)  K1 


y.  >  0 

7 1  — 


u  e  U  , 


0  V  x.  e  L. 
l  f 

V  x.  e  B. 

i  f-1 

V  x,  e  L, 

l  f 


where  L^Cfx^/i  e  N)  . 


If  in  problem  I  exists  a  basic  optimal  solution  with  {y^  =  0, 

V Xj  e  L^/L^},  then  all  the  basic  optimal  solutions  of  II  are  basic  optimal 
solutions  of  I.  Moreover,  for  the  networks  representing  the  above  problems, 
all  minimal  cuts  of  II  are  minimal  cuts  of  I. 


Proof 

Take  the  solution  of  I  with  (y^  =0,  Vx^  e  This  solution 

satisfies  the  constraints  of  II  and  clearly  Max  F  =  Max  F  Now  take  all 

u£U  1  ueU 

basic  optimal  solutions  of  II  and  add  new  variables  {y^  ^  0,  Vx.  s  L£/Lf } 
as  non-basic  ones  with  value  zero.  Clearly  all  these  solutions  satisfy  the 
constraints  of  I.  Moreover,  since  Max  =  Max  F^,  and  the  new  variables 
were  added  as  non-basic  ones,  the  solutions  are  also  basic  optimal  solutions 
of  I. 

For  the  corresponding  networks,  take  now  all  solutions  of  I  with 
{y_.  =0,  V  x.  eLj/L^}.  Removing  the  links  corresponding  to  (yVx^  e  LJ/L^} 
will  give  the  network  configuration  of  II,  and  hence  all  the  minimal  cuts  of 
II  are  minimal-cuts  of  1. 

□  Theorem  2.1. 

Theorem  2.2 

With  the  notation  of  theorem  2.1,  suppose  that  there  is  no  set  L£ 
such  that  L£,  and  that  all  basic  optimal  solutions  for  Lf  are  also 

basic  optimal  solutions  for  L£.  Then  a  minimal  cut  for  is  CX, X) ^  where 

X  »  { (ieW)U  {s)/x^  e  L^}  . 

Proof 

The  proof  will  be  carried  out  by  contradiction.  Clearly  every  node 

{i/x^ e  L^}  belongs  to  X  since  there  are  no  upper  bounds  on  the  variables 

y  .  Now,  suppose  that  there  exists  a  minimal  cut  (X',X')  such  that 
i 

X'  *  XUj  where  x^  t  L^.  Then  we  can  add  to  the  network  of  problem  II  a  new 
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non-basic  variable  y.  _>  0  with  value  zero  (i.e.,  a  link  connecting  the 
source  with  node  j  with  zero  flow).  By  Theorem  2.1,  all  basic  optimal 
solutions  of  II  are  also  basic  optimal  solutions  of  I,  where  L£  =  L^Ux^ 
in  contradiction  to  the  assumption  that  there  exists  no  such  set. 

□  Theorem  2 . 2 

Corollary  2.1 

Under  the  assumptions  of  Theorem  2.2,  when  maximal  flow  is  achieved, 
all  links: 

{  (i,k)  e  L/x.  £  Lf ,  xkj£Lf>  , 
have  a  flow  equal  to  c^  and  all  links: 

{ (j ,i)  e  L/x.  e  Lf,  x^  i  Lf}  , 

have  flow  with  value  zero.  (Follows  from  Theorem  2.2  and  the  Max-Flow-Min-Cut 
Theorem) . 

□  Corollary  2.1 

Now  we  are  able  to  present  a  method  for  solving  problem  (2.8).  Re¬ 
call  that  step  1  of  the  algorithm  calls  for  solving  a  series  of  linear  pro¬ 
grams:  First  the  linear  program  corresponding  to  all  states  leaving  the 

boundary,  then  those  corresponding  to  all  combinations  of  (n-1)  states 

1  2 

leaving  the  boundary,  etc.  In  other  words,  if  and  L^  are  two  sets  of 

states  leaving  the  boundary  corresponding  to  two  successive  substeps  of 

2  1 

step  1,  then  cardinality  L^  <_  cardinality  L^.  Now  suppose  that  the  current 
substep  of  step  1  corresponds  to  the  set  of  states  L^.  If  in  an  earlier 
substep  (say,  corresponding  to  the  set  Lp  we  obtained  a  basic  optimal  solu¬ 
tion  with  (y^  ■  0,  Vx.  e  LpLf}  and  Lfc  Lp  then  by  Theorem  2.1  all 
optimal  solutions  of  the  current  substep  (Lp'  are  included  in  the  set  of 
optimal  solutions  of  the  earlier  substep  (Lp .  Therefore  we  do  not  need  to 
solve  the  linear  program  corresponding  to  the  current  substep.  In  algorithmic 


terms,  there  are  two  ways  for  checking  the  existence  of  such  a  set  L£: 


1.  To  search  all  solutions  obtained  in  earlier  substeps  corresponding 
to  sets  where  L^CL^,  while  looking  for  a  solution  with 
(y.»0Vxj£  Lf^Lf } * 

2.  To  find  the  Maximal  Flow  corresponding  to  the  current  substep  (L^) 

and  then  to  look  for  earlier  substeps  such  that  L^C  L£  with 

the  same  value  of  Maximal  Flow. 

Method  (2)  follows  from  the  proof  of  Theorem  2.1. 

Suppose  now  that  there  is  no  set  L£  such  that  L£  and  such 
that  the  problem  corresponding  to  has  an  optimal  solution  with 

{y.  -  0  VXj  e  L^/L^}.  In  this  case  we  apply  to  the  network  corresponding  to 
Lf  the  Maximal  Flow  Algorithm.  According  to  Corollary  2.1,  when  maximal 
flow  is  achieved,  the  network  can  be  partitioned  into  two  subnetworks  as 
shown  in  Figure  2.4: 


Figure  2.4  -  The  two  subnetworks  separated  by  a  minimal  cut. 


Recall  that  since  we  are  building  a  feedback  control  region,  we  are 

searching  for  all  extremal  values  of  (y^/x^  6  and  the  corresponding 

controls  {u../(j,k)  eL>.  But  we  notice  that  only  pivots  on  the  subnetwork 
J  K 

corresponding  to  X  will  lead  to  new  extremal  values  of  {y^/x  6  so 


that  in  order  to  obtain  all  optimal  solutions  to  the  problem  corresponding 


to  we  concentrate  on  the  reduced  network  depicted  in  Figure  2.5: 


Figure  2.5  -  Network  configuration  before  the  pivoting  operation. 


In  Figure  2.5  we  define: 


'id 


1  < 

jeE(i) 

y  Lf 


V  is  1(d)  n 


b. 

l 


l  c.. 

J«H(i) 

V  Lf 


Cid 


Vifl(d)  and  i  e  I  (k/^  e  Bfl)  (2.10) 

Vic  1(d)  and  i  t  I  (k/^  e  Bfl) 


0 


otherwise 


I  b- 


(2.10a) 


xieLf 


For  the  first  basic  optimal  solution  we  take 


y.  *  b.  V  x.  e  L- 
'i  l  if 


u„  *  0  V  e  Lf,  j  e  E(i)  0  (i/x.  e  Lf} 


(2.11) 


By  pivoting  on  the  network  of  Figure  2.5,  we  obtain  all  the  required 
optimal  solutions,  when  for  each  solution  the  flow  values  on  the 'links  in  X 
are  equal  to  those  achieved  when  applying  the  Maximal  Flow  Algorithm. 
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Notice  that  at  every  substep  of  step  1  a  further  reduction  of  the  net¬ 
work  is  obtained.  Moreover,  from  considerations  of  Theorem  2.1,  there  might  be 
many  substeps  for  which  there  is  no  need  to  solve  a  linear  program. 


D.  States  Leaving  the  Boundary  from  Constructed  Feedback  Control  Regions 

We  deal  here  with  step  2  of  the  algorithm  of  Part  A  of  this  Section. 

In  this  step  we  find  the  leave- the- boundary  values  of  the  costates  for  all 

x^  e  8^  corresponding  to  a  constructed  feedback  control  region  R^,  and  then 

we  solve  a  linear  program  of  the  form  (2.2)  for  each  L  c.  B  .  Step  2  is  per- 

P  P 

formed  iteratively  for  each  previously  constructed  feedback  control  region.  We 
point  out  here,  that  the  values  of  the  costates  corresponding  to  I  and  the 
values  of  the  costates  corresponding  to  states  leaving  the  boundary  (L^)  are 
positive.  These  details  will  be  proved  in  Section  3.  As  for  step  1  of  the 
algorithm,  step  2  can  also  be  greatly  simplified,  as  shown  by  the  following 
theorem: 

Theorem  2.2 


s.t. 


Consider  the  problem 

Maximize  \  a.y.  a.  >  0  V  x.  e  I 

x.  el  .  1  1  1  i  P-1 

i  p-1 


'  l  uii+  £  \i  +  yi*°  VxieIp-l 
jeE(i)  13  kel(i)  1  1  p  x 

T  u.  .  -  l  u.  .  *  0  V  x.  e  B  , 

jeE(i)  13  kel(i)  1  P_1 

y.>°  V^elp.j 

ueU 


(2.12) 


(2.12a) 


A  basic  optimal  solution  of  (2.12)  with  arbitrary  coefficients  (a^  is  also 
a  basic  optimal  solution  of  (2.12),  with  (a^  »  l.Vx^el  j}. 
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Proof 

Assume  that: 

F.  =  optimal  value  of  (2.12)  with  {a.  =  1,  Vx.  el  .}  . 

1  l  l  p-i 

y*  »  optimal  solution  vector  to  (2.12)  with  (a.  >  0,  V  x.  e  I  ,} 
-  i  l  p-1 

Now  suppose  that: 


n  * 

l  y.  <  f, 

H  x  l 


(2.13) 


x.  el  , 
l  p-1 


Then  the  Maximal  Flow  Theorem  (see  Appendix  A)  implies  that  if  (2.13)  is 
satisfied, -we  can  always  find  a  path  between  the  source  and  the  destination 
nodes  of  the  network  corresponding  to  problem  (2.12)  on  which  we  can  increase 
the  flow,  in  contradiction  with  the  assumption  that  y*  is  an  optimal  solu¬ 
tion  to  (2.12)  with  (ai  >  0,  Vx^I^).  Now*  noting  that: 

l  X-  >  F  >  (2.14) 

x.el  , 
i  p-1 

cannot  be  satisfied,  since  clearly  a  flow  satisfying  (2.14)  violates  the  con¬ 
straints  (2.12a)  by  the  Maximal  Flow  Theorem,  then: 


[  y*  =  F  .  (2.15) 

x.el  ,  1  1 

l  p-1 

Moreover,  since  the  structure  of  the  network  corresponding  to  problem  (2.12) 
does  not  depend  on  the  weightings  {a^},  then  all  basic  optimal  solutions 
of  (2.12)  with  (a^  >0,  V  e  I  are  also  basic  optimal  solutions  of 

(2. 12)  with  (a.  *  1,  V  x.  e I  , }. 

i  i  p-i 

□  Theorem  2 . 2 

Now,  transforming  the  linear  program  (2.2)  into  a  Maximal  Flow  Prob¬ 


lem  in  a  network  with  weights  on  the  links  connecting  the  source  with  every 
node,  we  have: 
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s.t. 


Maximize 


x.e 

1 


I  uii  +  l  *  Yi  *  0 
jeE(i)  13  kel(i)  K1 

I  u. .  -  l  u,  .  =  0 
jeE(i)  13  kel(i)  kl 

y .  >0 

i  — 

u  e  U 


V 

V 

V 


x.  e  I 

l 


x.  e  B 
l 


x.  e  I 
i 


P-1 

P-1 


P-1 


h 


(2.16) 


(2.16a) 


According  to  Theorem  2.2,  all  solutions  of  problem  (2.16)  are  solu¬ 
tions  of  problem  (2.9)  with  »  I  ^ ,  therefore  every  linear  program  of 
step  2  of  the  algorithm  of  Part  A  of  this  section  reduces  to  the  following 
simple  problem: 


From  among  all  solutions  of  problem  (2.9)  with  Lf  ■  I  j,  (denoted 
by  y*)  choose  those  satisfying: 

Max  £  \.y. 

ye(y*}  x.el^  1  1 

Note  that  step  1  is  performed  for  all  possible  combinations  of  states,  assur¬ 
ing  that  we  can  always  find  a  set  such  that  *  I 
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Section  3 


THEORETICAL  RESULTS 


A.  The  Structure  of  the  y-Constraint-Figure 


in  the  Positive  Orthant  of  the  y-Space 


In  [M3]  the  y-constraint-Figure,  /,  is  defined  as  follows: 
V  -  (y  e  Rn/u  e  U  } 


where 


0  <  u.,  <  c.. 
—  lk  —  ik 


l‘/(i,h)  e  L 

and  the  linear  transformation  relating  Y  and  U  is  given  by: 


y(t)  =  -xCt)  *  -Bu(t)  , 

where  B  is  the  incidence  matrix  of  the  network. 

Since  U  is  a  bounded  convex  polyhedron  in  Rm,  its  image  /  is 
clearly  a  bounded  convex  polyhedron  in  Rn.  Thus,  every  face  of  Y  can  be 
analytically  described  by  an  expression  of  the  form: 


7  a.y.  =  i({a. })  , 

i,l  11  1 


(3.1) 


where  (a. >  is  a  set  of  coefficients  and  L({a. })  a  constant  with  value 
i  i 

depending  on  the  set  {a^,}. 

Our  first  aim  is  to  prove  that  all  coefficients  {a^}  corresponding 
to  any  face  of  Y  that  is  in  the  positive  orthant  of  the  y-space  are  nonnega¬ 
tive.  This  is  easily  demonstrated  if  we  consider  a  constraint  of  the  form 


n 

l  ay  <  L({a  }) 
i-1  1  1  1 


(3.2) 
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passing  through  the  positive  orthant  of  the  y-space  so  that  there  is  at  least 
a  point  y*  lying  on  the  boundary  of  (3.2),  satisfying: 

yj  >  0  Vi  e  M  .  (3.3) 

Suppose  now  that  a  coefficient  of  (3.2),  say  a^,  is  negative  and  let  us  con¬ 
centrate  on  one  of  the  paths  carrying  flow  from  the  node  k  to  the  destination. 
At  least  one  such  a  path  exists  since  y*  >  0  by  (3.3).  Clearly,  by  decreas¬ 
ing  the  flow  on  this  path  the  flow  on  the  network  will  remain  feasible,  but  we 
violate  the  constraint  (3.2)  since  a^  is  negative.  Therefore  all  coeffici¬ 
ents  of  (3.2)  must  be  nonnegative.  Moreover,  clearly  the  constant  Jl({a^}) 
is  given  by: 

n 

Ufa*})  =  Max  Y  a.y.  .  (3.4) 

y  i-i  1  1 

In  order  to  obtain  an  explicit  set  of  constraints  defining  V  in  the  positive 
orthant  of  the  y-space  we  must  take  all  the  constraints  (3.2)  and  eliminate 
those  which  are  redundant.  To  do  this,  we  shall  present  now  three  redundancy 
conditions  and  we  will  prove  that  if  a  constraint  such  as  (3.2)  satisfies  at 
least  one  of  the  conditions,  it  is  redundant.  Moreover  we  shall  prove  that 
a  constraint  satisfying  not  one  of  the  above  conditions  is  not  redundant.  In 
other  words,  a  necessary  and  sufficient  condition  for  a  constraint  of  the  form 
(3.2)  to  be  redundant  is  the  fulfilment  of  at  least  one  of  the  conditions. 

For  a  constraint  such  as  (3.2)  we  denote 
D  =  {i/0  <  a.  e  (a. )} 

l  l 

k(D)  »  Max  l  y. 

y  D  1 


and 
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Redundancy  Condition  1 
The  hyperplane 

l  a  y  =  i(ta.»  (3.5) 

D  1  1 

is  redundant  if  at  least  two  of  its  coefficients  are  different. 

Proof 

By  Theorem  2.3: 

{ye  V/l  a  y.  =  5.({a  })}  G  {ye  Y/\  y  =  k(D)}  . 

D  D 

Namely,  to  a  hyperplane  of  the  form  (3.5)  always  corresponds  another  hyper¬ 
plane  with  unity  coefficients  containing  all  points  of  tangency  of  (3.5)  with  Y. 


Redundancy  Condition  2 
The  hyperplane 

l  y.  »  k(D) 

D  1 

is  Tedundant  if  there  exists  a  set  of  indeces  D'  such  that: 


with 


D  C  D’ 

I  =  k(D) 
D'  1 


Proof 

By  Theorem  2.1,  if  (3.7)  is  satisfied,  then: 


(3.6) 


(3.7a) 

(3.7b) 


{y  e  Y/l  y .  ■  k (D) }  C  (y  e  K/T  y  =  k (D) }  . 

D  1  D' 

Namely,  the  hyperplane  (3.7b)  contains  all  points  of  tangency  of  (3.6)  with  Y. 


I  y.  -  k(D) 
D  1 


Redundancy  Condition  3 
The  hyperplane 


(3.8) 


is  reuundant  if  there  exist  sets  of  indeces  (B^,  j  =  1,2,  ...r},  satisfying 


l  k(B  )  =  k(D)  , 


(3.9a) 


J  =  1 
r 


U  B  -  D 
j=l  3 


>  (3,9b) 


5.  n  B.  =  1>)  Vj,,j0  e  (1,2,..  ,,r}  | 

h  J2  12  ; 


where 


Proof 


k(B  )  =  Max  l  y 
3  V  B.  1 
J 


(3.9c) 


There  are  no  common  links  to  the  minimal  cuts  corresponding  to  the 

hyperplanes  (3.9c),  because  the  existence  of  one  such  link  implies: 

r  r  r 

k(D)  *  Max  l  y.  =  Max  l  l  y.  <  [  Max  l  y  =  l  k(B,) 

V  D  1  -  ]=1  B.  1  j=l  y  B.  j=l  J 


J 


J 


or  k(D)  <  2  k(B .)  contradicting  (3.9a) 


J-l 


J 


Now  if  there  are  no  common  links  to  the  minimal  cuts  corresponding 
the  hyperplanes  (3.9c),  then  by  (3.9b)  follows  that: 

r 

{y  t  V/i  y.  =  k  CD) }  f  (ye //ye  O  (I  Y  ■  =  k(B  ))}  , 

-  D  '  ‘  3-1  »j  J 

therefore  (3.8)  is  redundant. 

□  Redundancy  Conditions 


We  prove  now  by  contradiction  that  a  hyperplane  of  the  form 

l  X:  =  k(D) 

rv 


(3-10) 


satisfying  none  of  the  three  redundancy  conditions  is  not  redundant.  To  this 
end  we  assume  that  (3.10)  does  not  satisfy  redundancy  conditions  2  and  3 


34 


(condition  1  is  clearly  not  satisfied  by  (3.10),  but  is  redundant.  If  (3.10) 
is  redundant,  then  there  is  a  non-redundant  hyperplaue  of  the  form: 


such  that: 


ly,  -  k(D') 
D'  1 


(3.11) 


(ye  y/l  V.  =  k(D)}  C  (ye  Y/  £  y.  =  k(D')>  , 
D  1  D’ 


(3.12) 


and  this  for  the  following  reasons: 


Every  solution  of  (3.10)  over  V  is  feasible  and  tangent  to  The 

surface  of  tangency  common  to  (3.10)  and  Y  is  clearly  a  convex  hyperplane  of- 
dimension  m,  where  m  <  (n-1),  (if  m  =  (h-1)  then  (3.10)  is  not  redundant). 
Now,  an  m-dimensional  hyperplane  on  the  boundary  of  V  corresponds  to  the  inter 
section  of  (n-m)  non-redundant  hyperplanes  of  dimension  (n-1),  and  every 
point  of  tangency  of  (3.10)  with  Y  belongs  to  every  one  of  these  (n-m) 
hyperplanes.  Hence,  (3.12)  applies  to  every  one  of  the  hyperplanes  forming  the 
surface  of  tangency,  and  for  the  proof  we  need  to  consider  only  one  of  them, 
say  the  hyperplane  given  by  (3.11). 

Now,  we  prove  that  (3.10)  is  not  redundant.  Note  that  in  general  the 
index  sets  D,  D'  can  be  partitioned  as  follows: 

D  =  (D1,D2> 

D'  =  (Dj.Dj) 

where 

Dl°  °2  “  V  °3  =  °2A  °3  =  {*}  • 

Our  aim  is  to  contradict  (3.12),  considering  all  possible  relations  between 
the  sets  D  and  O'.  There  are  four  cases: 
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Case  1:  Dj  *  {$};  D 2  +  {<(>};  ?  {<()}. 

In  this  case  we  choose  a  point  y*  satisfying: 

I  y?  =  =  UD)  , 

d  1  d2  1 

where 

y*  >  0  V  i  e  D_ 

J  l  2 

y*  =  0  Vi  e  D.  . 
i  3 

Clearly,  since  /  {$}  we  have: 

ly*  =  l  y?  =  0  <  k(D’)  , 

D'  °3 

thus  contradicting  (3.12). 


Case  2:  D3  =  {$);  Dx  /  {$};  D2  /  {*}. 

Here  we  choose  a  point  y*  satisfying: 

l  y|  =  k(D) 

D  x 


l  Yl  *  k(D  ) 
D„  1  * 


But  if  (3.10)  does  not  satisfy  redundancy  condition  3,  then 


k(D)  <  k(Dj)  ♦  k(D2)  , 

therefore  from  (3.13)  and  (3.14)  follows  that: 


l  y?  ■  lyf  •  k(D) -k(D,)<k(D.)  -  k(D') 
D'  1  Dl  1  *  1 

or 

<  k  (D  * )  , 

D’  1 


in  contradiction  to  (3.12)- 
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Case  3 

Let  us  choose  a  puj.nt  »’  satisfying  }  y'  =  k(DI  =  k(D  > 

6  1  1 

where 


y*  0  V  i  r.  o  *  D 
1  i  f  i 

y *  -  o  V  .  t  Dt 

A  ^ 

Now,  if  \  y  ’  =  k  CD  * )  -  k  (0 .  j ,  then,  situe  CCO',  hyperpiane  (3  10) 
D '  1  x 

satisfies  redundancy  condition  2,  contradicting  the  assumptions,  thus. 


I  y*  -  k(D')  , 

O' 


m  contradiction  to  (3.12). 

Case  4:  ?  { <p / ;  D2  t  1$};  f  l }  - 

In  this  case  we  choose  a  point  satisfying: 


where 


l  r  »  k(D)  i 

D  1  i 

I 

! 

v  *  >  0  V  i  ►  D 

i  i 

y*  *  0  VicDj  j 


This  point  satisfies  £  y*  -  k(D')  since  if 

D'  1 


(3  15) 


Iy?*kiD')  ,  (3.16) 

D'  A 

then  from  (3- IS)  and  (3.16)  follows  that: 


Max  l  y  -  l  y*  -  K (D ) 
V  DVD’  1  D'UD  1 


Thus,  (3.10)  satisfies  redundancy  condition  2  in  contradiction  to  the  assumptions 
Therefore  y*  contradicts  (3  12). 

% 

In  order  to  illustrate  the  procedure  for  obtaining  the  y -constraint 


figure  in  the  positive  orthant  cf  the  y-space,  we  present  the  fci rowing  .simple 


example  in  two  dimensions: 


Consider  the  network  depicted  in  Figure  3.1: 

c  -2 
12 


C21  *  1 


c  =2 
Id 


c  *  2 
2d 


Figure  3.1  -  Network  of  Example  3.1 


Finding  the  maximal  flow  values  of  the  networks  shown  in  Figure  3.2, 


we  obtain: 


(a)  Max  y.  =  4 

y  l 


(b)  Max  y2  =  3 


(c)  Max  (y.+yJ  *  4 

[/  1  4 


Figure  3.2  -  Networks  to  find  the  maximal  flow. 


CleaTly  the  hyperplane  yj  >  4  is  redundant  by  redundancy  condition  2,  and  V 


is  defined  by 


The  region  V  in  the  positive  orthant  of  the  y- space  is  depicted  in 


Figure  3.3: 


Figure  3.3  -  Structure  of  Y  in  Example  3.1 


□  Example  3.1 

B.  Special  Properties  of  the  Feedback  Solution  for  Single  Destination  Networks 
with  all  Unity  Weightings  in  the  Cost  Functional 

As  we  mentioned  in  Part  A  of  Section  2,  several  features  associated  with 
the  Constructive  Dynamic  Programming  Algorithm  in  the  case  of  multiple  destina¬ 
tion  networks,  or  when  the  cost  functional  has  nonequal  weightings,  complicate 
the  formulation  of  a  computational  scheme  to  implement  the  algorithm.  However, 
in  the  case  of  single  destination  networks  with  all  unity  weightings  in  the  cost 
functional  various  simplifications  result,  thus  permitting  us  to  develop  a  com¬ 
pact  algorithm  for  building  the  feedback  solution.  These  simplifications  are 
stated  and  proved  in  [Ml].  However,  the  approach  in  this  work  permits  us  to 
carry  out  the  proofs  in  a  more  simple  and  straightforward  manner. 

We  begin  with  the  central  theorem  of  this  section.  The  proof  of  this 
theorem  is  the  basis  of  all  subsequent  proofs. 

Theorem  3.1 

The  value  of  the  leave- the-boundary  costates  are  unique  at  a  given 
P  ‘ 


boundary  junction  time  t, 
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Proof 

Assume  that  t  is  a  fixed  boundary  junction  time  and  consider  first 

the  case  in  which  we  allow  all  states  lying  on  the  boundary  to  leave  it,  that 

is  L  *  B  .  Let  us  denote  by  t*  the  time  just  before  the  states  in  B 
P  P  P  P 

leave  the  boundary,  and  by  t  the  time  just  after  the  states  in  B^  have 
left  the  boundary,  (going  backwards  in  time).  Now,  the  restricted  Hamiltonian 


at  t  is: 
P 


H(t+) '•  z(t+)  *  T  A.  (t+)y.  , 

P  P  x  i  P 

i  P 


(3.17) 


and  the  set  of  operating  points  in  the  y-space  is  given  by: 


Y  =  {y  e  V7z(t*)  =  Max  £  A  (t*)y  >  . 

P  -  P  V  x.el  1  P  1 

i  P 


We  assume  in  this  proof  that  all  operating  points  in  are  nonnegative, 
that  is,  they  are  in  the  positive  orthant  of  the  y-space.  Later  we  shall  prove 
that  in  the  case  of  single  destination  networks  with  all  unity  weightings  in  the 
cost  functional,  we  can  consider  only  this  orthant  and  that  there  always  exists 
a  solution  there. 

We  denote  by  "a"  the  optimal  value  of  the  restricted  Hamiltonian  (3.17), 


that  is: 


Max  [  y.  » 
y  xieIp 


(3.18) 


By  the  Constructive  Dynamic  Programming  Algorithm,  we  have  to  find  the 
leave -the- boundary  costate  values  associated  with  all  states  in  B^  which  allow 
for  the  optimal  departure  of  the  states  from  the  boundary  at  tp,  (see  [M2]). 

In  geometrical  terms,  the  above  step  calls  for  finding  the  values  of  (A  (t  ), 

J  r 

Vx,  e  B  ),  for  which  the  global  Hamiltonian 
J  P 


H(t-):  z 


(t”)  *  I  A  (t+)y.  ♦  T  A.  (t")y. 
P  x.£l  1  P  1  x.iB.  1  P  1 


(3.19) 


i  P 


i  P 


contains  a  Bp-positive  face  of  /  (see  [M3]).  Notice  that  since  the  operation 
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is  made  by  rotating  H(t  )  around  H(t  ),  then  Y  e  H(t) , 

P  P  P  F 


Max  l  A  (t *)y.+  l  A  (t')y  «  Max  £  A(t*)y 

-TiPJ-vs-R1?1  V  •*  ,-T  1  P  1 


y  x.  el 
i  P 


x.  eB 
i  P 


Y  x.  el 
i  P 


Now  we  prove  the  assertion: 


therefore: 

*  a  .  (3.20) 


Assertion  3.1 


(i) 

A . (t+)  >  0 
l  P  ~ 

Vx.  e  I 
i  P 

(ii) 

A.  (t*)  >  0 
l  P  ~ 

Vx.  eB  . 
i  P 

Proof 

• 

(i) 

If  x.  e  I  , 
i  P 

then  x^  has  left  the  boundary  for  the  last  time  at  x. 

where  t  >  t  with  an  optimal  slope  >  0.  Now,  clearly  A^t)  >_  0 
since  otherwise  it  is  not  optimal  for  x ^  to  leave  the  boundary  at  t. 
Therefore,  since  ^(t)  *  -1,  Vte  (tp*T),  we  have  that  A^t*)  >  °» 

VxieIp’ 

(ii)  If  A. (t~)  <  0  for  some  x.  e  B  ,  then  it  is  not  optimal  for  x.  to 
l  p  I  P  1  . 

leave  the  boundary,  but  x^  does  leave  the  boundary,  so  that 
A.  (t")  >  0,  Vx.  e  B  . 

i  p  -  l  p 

□  Assertion  3.1. 

Now,  since  (A^t*)  >  0,  Vx^I^},  then  by  Theorem  2.3  we  have  that: 


(y  e  y/  l  A  (t*)y.  *  a}  C  {y e  Y/  j>  y.  *  k(I  )> 
x. el  1  p  1  *  x.el  p 


(3.21) 


i  P 


i  P 


wheT# 


k(I  )  *  Max  l  y  . 

P  1/  _  _  T  1 


*  V*p 


Moreover,  since  (A,(t‘)  >  0,  V  x.  eB  },  then  by  Theorem  2.3  it  follows  that: 
ip—  l  P 


-  41  - 


(y *17  I  ^(t*)y.  ♦  I  \(Oy.  *a>C{ye//  l  y .  ♦  [  y.  -k(l  UA)}, 

'  x.el  1  P  1  x.eB  1  p  1  ‘  x.el  1  x.eA  1  p 

1  p  X  P  X  P  1 

(3.22) 

where  A  is  a  set  of  states  such  that  AC  B  .  Clearly,  since  Y  eH(t"), 

P  P  P 

then  by  (3.21)  and  (3.22),  Y^  belongs  to  all  sets  in  (3.21)  and  (3.22), 
therefore : 

Max{  l  y.  ♦  l  y)  *  Max  l  y.  «  k(I  )  .  (3.23) 

y  x.el  1  x.eA  1  /x.el  1  p 

i  p  i  i  p 

Next,  we  prove  the  following  assertion: 

Assertion  3.2 

There  is  a  (unique)  maximal  set  A  satisfying  (3.23). 

Proof: 

Suppose  there  exist  m  sets  (A^ ,  j  »  1,2, . . . ,m),  satisfying: 

Max*  l  ♦  l  y.)  «  k(I  ),  j  *  l,2,...,m  . 

V  x.elp  1  x.eA.  1  p 

Now,  consider  the  following  networks: 

(a)  Network  having  links  connecting  the  source  with  every  node  in  1^, 

namely,  network  corresponding  to  the  problem  Max  £  y. ,  (see  parts 

y  x.el  1 

B  and  C  of  Section  2) ,  1  p 

(b)  Network  having  links  connecting  the  source  with  every  node  in  1^  and 

in  A^,  namely,  network  corresponding  to  the  problem 

Max  (  l  y.  ♦  I  y.). 
y  V!/  x.eAl 

(c)  As  (b),  but  Ax  being  replaced  by  ' 

Now  by  Theorem  2.1,  all  networks  in  (a),  (b)  and  (c)  have  at  least  one 
common  minimal  cut.  But  this  minimal  cut  is  obviously  a  minimal  cut  of  the 
network  corresponding  to  the  problem: 


so  that 


-  42  - 


Max(  l  y .  *  l  Yi ) 

V  x.el  1  ■  1 

p  x.  g  U  A 


1  j»l  3 


Max{  l  y  *  l  y.  >  «  k(I  ) 

1/  „  _T  1  ~  1  P 


V  x.el 
i  P 


x  e  U  A. 
1  j=l  3 


Therefore,  since  there  is  a  finite  number  of  states,  then  there  is  a  (unique) 


maximal  set  A  *  \j  A,  satisfying  (3.23), 


□  Assertion  3.2 


From  now  on,  the  face  of  V  (of  dimension  less  ot  equal  to  n-1), 

defined  by  the  set  I  U  A  will  be  referred  to  as  the  (I  UA)-face. 

P  P 

Now,  we  concentrate  on  the  set  of  states  Bp  and  we  prove  the  follow¬ 
ing  assertion: 

Assertion  3.3 

If  x  eB  but  x.  £  A,  then  X.(t")  *  °* 
j  P  1  IP 

Proof: 

The  proof  will  be  carried  by  contradiction.  Assume  that  (tp)  =  0 

for  a  set  of  states  D  such  that  DCB  /A.  Then,  the  enlarged  Hamiltonian 

P 

is  given  by: 

h<9:  i(9  ■ ,  l  JoVV'j  (3-24) 

i  p  i  j 

Moreover,  from* Theorem  2.3  follows  that: 

{yty/2(0  -  a)  C  (ye//  [  y  ♦  I  y.  ♦  I  y.  -  k(I  UAUD)}  .  (3.25) 

p  "  x.el  1  x.eA  1  x.eD  3  p 

i  p  i  j 

Now,  if  k(I  UAUD)  >  k(I  ),  then  Y  i  H(t*)  which  is  a  contradiction;  and 
P  P  P  P 

if  k (I  UAUD)  ■  k (I  ),  then  clearly  DCA  since  A  is  the  maximal  set 
P  P 
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satisfying  (3.23),  again  a  contradiction.  Therefore  (A.(t~)  *  0,  ¥  x.  e  B  /A) 

3  P  3  P 

□  Assertion  3.3. 


From  assertion  3.3  follows  that  the  enlarged  Hamiltonian  is  given  by: 


H<v:  j(v  *  x  L  ■  •  ■  (3-26) 

ip  i 


and  that  the  Bp-positive-face  is  given  by  the  intersection  of  the  Hamiltonian 
(3.26)  with  the  (IpU  A)-face  (3.23). 

Now  we  prove  by  contradiction,  that  the  Hamiltonian  (3.26)  is  unique, 
implying  that  the  Bp-positive-face  is  unique.  Suppose  that  there  are  two 
Hamiltonians  of  the  form  (3.26)  such  that  the  intersection  of  every  one  of 
them  with  the  (I  U  A) -face  gives  two  distinct  B  -positive- faces.  These  Hamil- 
tonians  are: 

z^t")  =  l  A.  (t+)y.  ♦  l  Al(t')y.  *  a  ,  (3.27a) 

P  x.  el  1  P  1  x  eA1  P  1 

l  p  l 

z2(t")  *  l  A. (t+)y.  +  I  A?(t")y.  *a  ,  (3.27b) 

P  x. el  1  P  1  x.cA  1  P  1 

l  p  l 


where  for  the  two  Hamiltonians  to  be  distinct,  there  is  at  least  one  x.  eA 
1  -  ? 

such  that  A^(t  )  /  Ai(t  ).  This  also  implies  that  there  exist  at  least  two 
1  2 

points  y  >  0  and  y  >  0  of  the  (IpUA)-face,  such  that: 


and 


zl(tp)  *  a 
z2(t‘)  <  a 
z2(t‘)  »  a 
z^t'}  <  a 


f  at 


(3.28a) 


(3.28b) 


Now  consider  the  intersection  of  one  of  the  Hamiltonians,  say  (3.28a) 


with  the  (IpUA)-face.  Since  both,  the  Hamiltonian  and  the  (IpUA)-face  are 
convex  and  linear,  their  intersection  will  give  a  convex  polytope.  Moreover, 
since  the  (I^U  A) -face  is  bounded,  the  intersection  will  give  a  convex  poly¬ 
hedron.  Thus,  starting  at  an  extremal  point  of  the  polyhedron,  we  can  always 
reach  any  other  extremal  point  of  the  polyhedron  through  a  series  of  pivots. 
Note  that  all  extremal  points  of  Y^  are  contained  in  the  above  intersection. 
Suppose  that,  starting  at  some  extremal  point  of  Y  ,  we  want  to  reach  the 
point  y1  through  a  series  of  pivots  in  such  a  way  that  after  every  pivot  of 
the  above  series  we  reach  an  extremal  point  of  the  polyhedron,  that  is,  satis¬ 
fying  the  equations  of  the  Hamiltonian  and  the  (I  U  A) -face.  We  state  that 

P 

the  transition  between  Y^  and  yL  can  be  done  in  the  following  way: 


(i)  Starting  at  Y^,  reach  the  extremal  point  y  through  a  series  of 
pivots  on  the  polyhedron,  where 

a 


ya  *  y1  V  x.  e  I  U  A 

1  1  1  P 


ya  =  0  V  x.  e  B  /A 

l  p 


(3.29) 


(ii) 


Without  changing  the  coordinates  corresponding  to  the  states  in 
reach  y1 . 


IpU  A, 


The  above  procedure  is  justified  by  noting  that  if  we  are  moving  over 
the  (IpUA)-face,  then  every  point  reached  after  a  pivot  operation  must 
satisfy: 

I  y •  ♦  l  yA  -  k(i  )  ,  (3.30) 

x. el  1  x.eA  1  P 

i  p  i 


and  by  Corollary  2.1  the  network  looks  as  in  Figure  3.4: 


on» 
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in  I  U  A  and  B  /A  because  by  making  such  a  pivot  we  violate  (3.30). 

P  p 

Recall  that  our  aim  is  to  contradict  (3.27),  that  is  to  prove  that 
the  Hamiltonian  is  unique.  To  this  end  we  shall  deal  with  pivots  in  the  part 
of  the  network  corresponding  to  states  in  I^U  A.  First  we  prove  the  follow¬ 
ing  assertion: 


Assertion  3.4 

Starting  at  Y^,  we  can  always  find  a  path  in  the  network  between 
each  node  {j/x^  c  A}  and  some  node  {i/x.^  e  1^},  that  permits  us  to  perform 
a  pivot  while  increasing  (from  zero)  the  flow  in  the  link  corresponding  to  y 
and  still  to  remain  on  the  Hamiltonian  (3.27a). 


Proof 


Consider  the  series  of  pivots  that  reaches  y  from  Y  .  Suppose  we 


show  that  after  completing  n  pivots  of  this  series,  we  could  also  have  in¬ 


creased  the  flow  y^ ,  where  x.  eA,  corresponding  to  the  (n+l)-st  pivot 
before  the  n-th  pivot,  reaching  in  this  way  a  point  yb  that  satisfies  the 


equation  of  the  Hamiltonian  (3.27a).  Clearly,  the  application  of  the  same 
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arguments  to  the  new  series  of  pivots  that  reaches  y*5 
us  with  the  induction  step  that  proves  our  assertion. 


from  Yp  provides 


y(J) 

yl 

yfi) 

y2 

x<« 

y(n) 

Ay 


For  the  sake  of  clarity  the  following  notation  will  be  used: 

=  the  slack  variable  y.  that  increases  (from  zero)  in  the  i-th  pivot. 

»  • 

=  the  slack  variable  that  decreases  its  flow  in  the  i-th  pivot. 

=  the  costates  corresponding  to  the  slack  variables  of  the  i-th  pivot. 

=  the  point  reached  after  the  n-th  pivot. 

=  the  increase  of  flow  in  a  slack  link  after  pivoting. 


Suppose  we  cannot  perform  the  (n+l)-st  pivot  before  the  n-th  one. 
Clearly,  this  occurs  only  when  there  is  at  least  one  link  common  to  the  n-th 
and  (n+l)-st  pivot  cycles,  such  that  the  execution  of  the  (n+l)-st  pivot  is 
possible  only  after  the  n-th  one.  All  such  "critical"  links  must  be  in 
opposite  directions  on  the  n-th  and  the  (n«-l)-th  pivot  cycles.  We  shall 
contradict  our  assumption  by  showing  that: 

(a)  .  A(n)  . 


(b)  before  the  n-th  pivot  there  exists  a  cycle  that  permits  increasing 
y^n+1^  while  decreasing  . 

Consider  first  the  situation  before  the  n-th  pivot,  (see  Figure  3.5). 
Suppose  that  (k,j)  is  the  first  critical  link  on  the  (n+l)-st  pivot  cycle. 
Clearly,  since  the  link  is  critical  it  belongs  also  to  the  n-th  pivot  cycle. 
Then  follow  the  (n+l)-st  pivot  cycle  from  y^n+1^  up  to  the  node  k  and 
then  the  n-th  pivot  cycle  up  to  y ^ .  Clearly  we  can  perform  a  pivot  on 
this  cycle.  Moreover,  X^n+^  <_  X^  because  if  X^n+^  >  X^  then,  per¬ 
forming  the  pivot  corresponding  to  this  cycle  will  give: 
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l  A(t)y(n-1)  +  l  X  (t  )y. (n-1)  +  A(n  1  Ay  +  A  n^Ay  =  a  ♦  (X  n+1^  -  A  n  )Ay  ■*  a, 

x.el  1  p  1  x.eA  1  p  1 

l  p  x 

thus,  contradicting  the  optimality  of  the  Hamiltonian  (3.27a). 

Next,  consider  the  situation  after  the  n-th  pivot,  (see  Figure  3.5). 

Let  (k’.j1)  be  the  last  critical  link  on  the  (n+l)-st  pivot  cycle.  Consider 
the  following  cycle:  start  at  ,  follow  the  n-th  pivot  cycle  backwards 

up  to  the  node  k'  and  from  there  on  follow  the  (n+l)-st  pivot  cycle  (for¬ 
wards),  up  to  y^n+1^ .  Clearly  we  can  perform  a  pivot  on  this  cycle  and  hence 
A^  <_  A^n+1^  for  the  same  reasons  as  before. 

- n-th  pivot  cycle  links;  -  -  -  (n+l)-st  pivot  cycle  links;  ...  critical  links 

~\ 


Figure  3.5  -  Example  of  critical  links. 

Therefore  we  conclude,  that  A^  =  A^n+^  and  as  said  already,  before 
the  n-th  pivot,  there  is  a  cycle  including  y|n+1^  and  y^  on  which  a  pivot 
can  be  performed  so  that  we  can  increase  y^  J  (from  zero)  before  the  n-th 
pivot  while  reaching  a  point  that  satisfies  the  equation  of  the  Hamiltonian 
(3.27a).  This  proves  the  induction  step. 


□  Assertion  3.4 
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Now  we  contradict  (3.27).  By  assertion  3.4,  starting  at  Y  ,  we  can 
reach  points  y*  by  performing  a  single  pivot  on  the  network.  These  points 
lie  on  the  (I^U  A)-face  and  also  satisfy  the  Hamiltonian  equation  (3.27a), 
namely: 


I  X.  (t*)y?  +  X1 (t~)y*  *  a,  Vx  eA  . 
x  el  1  P  1  m  P  m  m 

i  P 


(3.31) 


We  check  now  the  second  Hamiltonian  (3.27b)  at  the  above  points. 

Notice  that  since  we  assumed  that  (3.27b)  is  optimal  over  Y,  it  must  satisfy: 

(3.32) 


F  X.  (t+)y7  +  X  it  )y*  <  a,  Vx  eA  . 
x  Jj  x'-  p^i  m  p  7m  -  m 

i  P 


From  (3.31)  and  (3.32)  follows  that: 


X*(0  <  xV),  Vx  eA  . 

m  p  —  m  p  m 


(3.33) 


But  obviously,  assertion  3.4  holds  also  for  the  second  Hamiltonian  (3.27b),  hence: 

(3.34) 


x!(tj  <  XV),  Vx  eA  . 
m  p  —  m  p  m 


Thus,  from  (3.33)  and  (3.34),  we  have: 

X*(0  =  eA  . 

m  p  m  p  m 

We  then  conclude  that  the  Bp-positive-face  is  unique  and  therefore  the  leave- 
the-boundary  costates  are  unique. 

When  only  some  of  the  states  leave  the  boundary,  we  consider  the  re¬ 
stricted  Hamiltonian: 


H  ,(t"):  z(t”)  *  y  X. (t+)y.  +  y  X.(t")y.  , 

P-1  P  P  x.e!  i  P  i  v  tr  1  P  1 


i  P 


x.eL 
i  P 


with  L  C  B  .  But  in  order  to  satisfy  the  necessary  conditions  the  optimiza- 
P  P 


tion  must  be  global,  namely: 

H  t  j  *  H(t“)n  * P  P  » 

P-1  P  P 


a  +p 


whete 
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a  *  cardinality  of  I 
P  P 

Pp  ^  cardinality  of  Lp  , 

V0D 

and  the  basis  vectors  of  1R  p  p  are  the  elements  of  Ip  Therefore  the 

values  of  the  leave- the-boundary  costates  in  the  case  of  L  £  B  are  identi 

P  P 

cal  to  those  corresponding  to  the  case  of  B  so  that  they  are  unique. 


□  Theorem  3.1. 


Theorem  3.2 

If  the  values  of  the  costates  corresponding  to  states  travelling  on 
boundary  arcs  during  the  interval  between  two  successive  junction  times 
tp+j  and  t  are  equal  to  the  leave-the-boundary  costate  values,  then 
every  costate  satisfies  exactly  one  of  the  following  conditions: 


A.  (t)  =  A.  (r)  =  0, 

l  l 

VTe(tp'tp+l)  ' 

(3.35a) 

A^t)  *  -1;  Ai(r)  >  0, 

Vte  (-,tp+i)  . 

(3.35b) 

Moreover,  in  every  feedback  control  region  there  always  exists  an  optimal 
control  satisfying  y  j>  0. 


Proof 

Every  optimal  trajectory  in  the  state  space  is  characterized  by  a 

control  switch  time  sequence  (t-.t-  t  ,,t  ,t  when  each 

r  r-i  p+i  p  p-i 

switch  time  represents  a  transition  between  two  neighboring  feedback 
control  regions  in  the  state  space  (see  [Ml]).  The  interval  between  two 
successive  switch  times  t  ,  and  t  is  characterized  by  the  set  of 

p+l  p 

states  travelling  on  interior  arcs,  Ip.  Considering  every  possible  control 
switch  time  sequence  insures  that  we  have  traversed  every  feedback  control 
region,  thus  covering  the  whole  state  space.  To  prove  the  Theorem  we  shall 
consider  a  typical  control  switch  time  sequence,  and  we  will  prove  by  induc¬ 
tion  that  the  costate  trajectories  satisfy  (3.35),  and  also  that  there  exists 
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at  least  one  optimal  control  foT  which  y  ^  0  within  every  feedback  control 
region.  Notice  that  the  existence  of  such  an  optimal  control  allows  us  to 
consider  only  the  set  of  trajectories  for  which  there  is  no  return  to  the 
boundary  of  states  (x  =  0),  backwards  in  time. 


Now  we  proceed  with  the  proof.  Consider  first  the  interval  (t  ,t^  ^) . 
Since  (X^(t£)  =  0,  Vx^cL^}  (see  [M2],  pp.  2S-26),  at  t£  all  costates 
corresponding  to  the  states  in  are  identical.  Now,  by  Theorem  3.1,  for 

the  costates  corresponding  to  states  on  the  boundary  to  be  at  their  leave- the- 
boundary  values,  we  have  to  find  the  maximal  set  of  states  Af  ^  that  satisfies 


Max{  I  y/  I  y  }  =  Max  £  y  =  k(L  )  .  (3.36) 

y  x.cL_  x.eA.  y  x.eL_  a 

if  l  f  - 1  if 


Moreover,  the  enlarged  Hamiltonian  at  t£  must  lie  on  /  in  such  a  way  that 
all  operating  points  in  satisfy  equation  (3.36).  But  as  said  before, 

all  costates  corresponding  to  states  in  l,£  have  the  same  value  at  tf,  so 
that  the  enlarged  Hamiltonian  will  contain  the  (l£  A^^-face  (3.36). 
Hence  Theorem  3.1  implies  that 


(i)  A^ft^)  =  0  for  a11  xi  e  Bf  l^f-i 


(ii)  the  values  of  all  X^(t^)  are  identical  for  all  i  such  that 

VLfUAf-i- 

where  (i)  follows  from  Assertion  3.3,  and  (ii)  follows  from  Assertion  3.4  and 
the  fact  that  all  costates  corresponding  to  states  in  [,£  have  the  same  value. 
Now,  clearly  there  exists  at  least  a  point  on  the  (1^  A^^-face,  satisfying 


(3.37) 


Moreover,  since  we  require  that  all  costates  corresponding  to  states  travel¬ 
ling  on  boundary  arcs  in  the  interval  (t^t^)  be  at  their  leave-the-boundary 


yi 

V 

o 

*• 

Vx.  e  L 

l 

yi 

*  0, 

Vxi  e  B 
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values,  (i)  and  (ii)  hold  for  all  xe  s0  that  the  point  (3.37)  is 

an  optimal  solution  for  all  xe  (t^t^  j).  Therefore  we  have: 

A.(x)  *  -1,  VxieLfUAf_1  =  If_1UAf  l 

^(x)  =  Ai(x)  =  0,  V  xi  e  Bf.1/Af_1 
Vxe  (tf,tf  l) 


(3.38) 


Now  we  perform  the  induction  step:  Consider  the  time  interval 
) ,  and  assume  that 

? 

least  one  point  satisfying: 


(t  +1»tp)»  and  assume  that  the  set  of  operating  points  Y^,  contains  at 


0,  V  x.  e  I 


y.  *  0,  ^x.  e  B 
i  l  p 

We  denote  by  Ap  the  maximal  set  of  states  satisfying 


(3.39) 


Max(  [  y,*  l  y. }  *  k(I  ) 


y  x.  el 
i  P 


x.  eA 
i  P 


Keeping  the  costates  corresponding  to  states  in  at  their  leave-the- 

boundary  values  at  every  instant  on  the  interval  will  result  by  Theorem  3.1  in: 

(3.40a) 


X.  (x)  *0,  Yx.  e  B  /A 
l  l  p  p 


Xi(x)  *  X^(x)  VxieAp,  xj  6  Xp 


(3.40b) 


Moreover,  from  the  existence  of  a  point  satisfying  (3.39)  during  the  entire 
interval  follows  that: 

■  -1.  Y*,el. 


Vte(VVi)  ■ 

therefore  from  (3.40)  and  (3.41)  we  have: 


(3.41) 


A 
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V(t)  *  X^t)  =  0,  )/^i  t  Bp/Ap  , 


(3.42a) 


^(t)  *  -1;  X^t)  >  0,  Vx^IpUAp  , 


(3.42b) 


V"('pV  - 

where  X^(x)  >  0  follows  from  Assertion  3.1  of  Theorem  3.1. 

Now  consider  the  boundary  junction  time  tp  where  we  allow  the  set 
of  states  Lp  to  leave  the  boundary,  so  that  1^  ;  s  IpU  Lp.  Noting  that 
we  have  kept  all  costates  corresponding  to  states  in  Bp  at  their  leave- the- 
boundary  values  during  all  the  interval  (t  ,t  ,),  we  realize  that  at  t 

p  p+1  p 

(3.40)  is  satisfied,  namely  at  the  junction  time  the  Hamiltonian  does  not 


change  its  orientation.  Therefore  in  order  to  reach  points  in  Yp  ^  (from 

a  point  in  Yp) ,  we  have  to  perform  pivots  on  the  subnetwork  corresponding  to 

states  in  I  U  A  (see  Figure  3.6),  increasing  in  this  way  the  flows  on  the 
P  P 

links  corresponding  to  y^,  where  x^  e  1^0  Ap,  and  then  to  increase  the  flows 

on  the  links  corresponding  to  y. ,  where  x.  eL  A  (B  /A  ). 

x  l  p  p  p 


u  ■  0 
ki 


B  /A 
P  P 


u  ■  c 
ik  ik 


i.6  -  Network  representing  transition  between  Y  to  Y  , 

—  p  p-1 
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As  a  result  of  the  increase  of  flow  on  a  link  y corresponding  to  an 

x^  e  Lpfl  Ap,  there  is  a  decrease  of  flow  on  some  link  y^  corresponding  to  an 

x.  e  I  ,  however,  it  is  clear  that  we  can  always  increase  the  flow  on  y.  in 
IP  i 

such  a  small  amount  that  y.  will  still  remain  nonnegative  for  all  x.  e  I  . 

J  IP 

Moreover,  the  increase  of  flow  on  links  y^  corresponding  to  x^  e  CB^/A^)0^p 

does  not  result  in  decreases  of  flows  on  other  slack  links  so  that  Y  ,  will 

P-1 

contain  points  satisfying 


y.  >  0  V x.  e  I  , 
x  —  l  p-1 


y.  »  0  Vx-  e  B  , 

v  l  p-1 


(3.43) 


Now  I  c I  ,  so  that  when  maximal  flow  is  achieved  on  the  network 
P  P-1 

corresponding  to  1^  ^  (i.e.,  network  with  slack  links  y^  for  all  x^  e  Ip  j 

only),  all  minimal  cuts  of  I  will  be  filled.  In  addition,  the  sets  I 

P  P 

and  I  U  A  have  the  same  minimal  cuts,  therefore: 

P 


(3.44) 


APC  Viulp  • 

where  Ap  ^  is  the  maximal  set  of  states  satisfying 

Max{  l  y.  +  l  y. )  *  k(I  )  . 

y  x.  el  .  1  x.  eA  .  1  P’1 

i  p-1  i  p-1 


Now,  from  the  existence  of  a  point  satisfying  (3.43),  and  from  the 

same  considerations  as  in  the  interval  (t  ,t  ,),  we  have 

P  P+1 


V'5  .  A.(t)  -  0  Vv'p-A-l 

^(t)  •  -1;  ^(t)  >  0;  Vxt  e  !p_iU  Ap  l 

V  r  e  (t  ,t  ) 

P-l  P 


(3.45) 


Moreover,  since  I  Cl  from  (3.45)  follows  that 
P  P-1 


XjCt)  ■  -1;  X^(t)  >  0;  Vxi  e  Ip  , 

V  t  e  (t  )  , 

P“l  P 


(3.46) 
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and  since  L  £l  , »  from  (3.44)  and  (3.45)  follows  that: 
P  P-l 

^CO  *  -1;  ^(t)  >  0,  VxieAp, 


Vxe  (Vi*V 

Therefore,  equations  (3 . 39) ,  (3.42),  (3.43), 
induction  step  that  proves  the  theorem. 


(3.47) 

(3.46)  and  (3.47)  provide  the 


□  Theorem  3.2. 


The  costate  trajectories  (3.35)  satisfy  the  necessary  conditions  (which 
are  also  sufficient)  derived  in  [M2].  Moreover,  these  trajectories  were 
obtained  by  checking  at  every  moment  the  global  Hamiltonian,  thus  insuring 
that  values  of  the  costates  corresponding  to  states  on  the  boundary  can  always 
be  found  such  that  the  solutions  obtained  by  considering  only  the  restricted 
Hamiltonian  are  also  globally  optimizing  solutions.  The  above  arguments  yield 
to  the  following  corollary: 

Corollary  3.2 

Any  solution  to  the  constrained  optimization  problem  (see  [M2],  p.  64) 
is  also  a  globally  optimizing  solution. 

□  Corollary  3.2. 


Theorem  3.3 

The  set  of  optimal  controls  does  not  switch  between  boundary  junction 
times;  that  is,  there  are  no  break  points  between  boundary  junctions. 

Proof 

Consider  the  time  interval  (t  ,t  ,)  and  denote: 

P  P+l 

Max  l  A  (t)y  -  a(t)  . 

v  x  cl  1  1 

r  l  p 

Now  consider  the  time  x.  e  (t  ,t  ,).  The  set  of  operating  points  is  given  by 

1  p  p*i 
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Y  CO  »  (y e  £  Ur.)y.  «  aCiJlHR 
p  1  “x.el  1  1  1  1 

i  P 
o_ 

where  the  basis  vectors  of  R  p  are  the  elements  of  I  . 

P 

Recall  that  all  points  in  (3.48)  satisfy 


P 


(3.48) 


Now  consider  the  time 
ing  points  is  given  by 


i  yi  *  k(y  • 

x.el  v 

i  P 


vyyp*  where 


(3.49) 

t^  +  At.  The  set  of  operat- 


a 

Y  (O  *  ty  e  V/  l  \  (t  )y.  *  a(x  ))n»p  .  (3.50) 

p  1  x.el  1  1  1  1 

i  P 

But  every  point  in  (3.50)  satisfies  (3.49),  and  also  by  equation  (3.42b)  we 
have 


A.(r)  *  -1,  Vxi  e  Ip 
V™  (tp,(tp+1)  . 

Therefore: 

Max  l  *4(Oyt  55  Max{  £  X .(rjy.  ♦  At  £  y . }  «  Max  £  X  (t  )y  *  At  k(I  )  . 
Y  x.el  121  y  x.el  1  1  1  x.el  1  V  x.el  111  P 

i  p  ip  ip  i  p 

Thus : 


W 


yv 


so  that  the  set  of  optimal  controls  does  not  switch  between  boundary  junction 
times . 

□  Theorem  3.3. 


Theorem  3":  4 

There  is  one  subregion  per  region  with  respect  to  any  set  of  state  vari 
ables  leaving  the  boundary. 

Proof 

Consider  the  feedback  control  region  R_  constructed  from  the  set  of 
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optimal  solutions  Y  on  (-*»,t  .).  Also,  let  L  be  the  set  of  state  vari- 

p  p+1  p 

ables  which  we  choose  to  allow  to  leave  from  R  backwards  in  time,  where 

P 

LpCBp.  According  to  the  definition  of  subregions  (see  [M2],  p.  57),  we  must 

show  that  the  state  variables  in  leave  with  the  same  set  of  optimal  controls 

from  every  point  of  Rp.  Recall  that  by  Theorem  3.3  the  optimal  trajectories 

associated  with  the  state  variables  in  Lp,  leaving  Rp>  do  not  experience 

a  break  as  time  runs  to  minus  infinity,  so  that  to  prove  that  there  is  one  sub- 

region  per  region  is  equivalent  to  showing  that  the  set  of  optimal  controls 

associated  with  the  state  variables  in  L  leaving  the  boundary  is  the  same 

P 

for  every  boundary  junction  time  t  e  (-*»,tp+j).  Note  that  we  need  not  be  con¬ 
cerned  with  the  common  breakwall  sequences  (see  [M2],  p.  56),  since  by  Theorem  3.3 
there  are  no  breakwalls.  Now  consider  two  different  boundary  junction  times. 


t  ,  t  e  (-«,t  .)  and  assume  that  t  >  t 


P2 


>1 


P2 


By  Theorem  3.3  we  have  that 


wv>  ■  {y-cV/.  £  Wyi  ’  a(t°-)in  * 


p-i'  pj 


Vp 


VlU  A 


(3.51) 


p-1 


and 


where 


a  +P. 


i  (t  )  ■  (y  e  *7  l  X  (t  )y  *  a(t  )}  f)  R  p  P 

“*  P-s  "  r  .  i  a  *  y<2  ro 


>-i^V 


p-i 


?2 


(3.52) 


a(x)  ■  Max 


/  I 


-I^Vl 


X. (f)y. 
l 


<3*0 


The  coordinates  of  R  p  P  are  the  elements  of  I 
as  in  Theorem  3.2. 


and  A 


p-1 


is  defined 


Now,  by  Theorems  3.2  and  3.3  we  know  that 

"  "l»  VxicIpUAp 
X^t)  -  0,  VxicBp/Ap 

VTe  (-.tp+1) 

where  A  is  defined  as  in  Theorem  3.2. 

P 


Hence: 


A.(t  )  =  A.(t  )  *  0,  Vx.  e  (I  ,U  A  )/(I  !JA  )  , 

x  p2  1  px  i  p-1  p-1  p  p 

*.(t  )  =  A  (t  )  +  (t  -t  ),  Vx.  c  I  UA  , 

1  p2  1  pl  P1  p2  p  p 


therefore 


Max  l  A  (t  )y  *  Max{  £  A  (t  )y  -  (t  -t  )  £ 

/  I  ,|JA  1  p2  1  y  I  ,1/A  1  P1  1  P1  p2  I  ,  Ij  A  . 

p-lw  p-1  p-ly  p-1  P-1  P-i 


Max  l  A  (t  )y.  +  (t  -t  )k(l  ) 

i  ViuVi  1  1  2 


so  that 


Y  (t  )  *  Y  (t  )  , 

p-1  px  p-1  p2 


then  the  set  of  optimal  controls  associated  with  state  variables  in  L 
leaving  the  boundary  from  Rp,  is  the  same  for  every  t  e  (-°°,t  ^)- 


□  Theorem  3 . 4 . 
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Section  4 


THE  ALGORITHM 


This  section  is  devoted  to  the  description  of  a  simplified  algorithm 
for  building  the  feedback  solution  to  the  minimun  delay  dynamic  message  rout¬ 
ing  problem  for  single  destination  networks  with  all  unity  weightings  in  the 
cost  functional.  The  simplifications  obtained  arise  from  the  special  proper¬ 
ties  of  these  kind  of  networks,  as  discussed  in  Section  3,  and  from  the  new 
approach  to  solve  the  linear  programs  presented  in  Section  2. 

In  Part  A  we  state  the  algorithm  while  explaining  its  steps.  In 
Part  B  we  present  a  method  for  obtaining  all  optimal  solutions  of  the  linear 
programming  problems  required  by  the  algorithm  and  finally  Part  C  brings  an 
example  that  provides  a  good  insight  of  the  performance  of  the  algorithm. 

A.  Statement  of  the  Algorithm 
Operation  1 

List  all  possible  trajectories  in  the  state  space,  by  writing  all  pos¬ 
sible  sequences  of  states  leaving  the  boundary  backwards  in  time. 

□  Operation  1. 

Recall  that  every  optimal  trajectory  in  the  state  space  is  character¬ 
ized  by  a  control  switch  time  sequence.  Each  switch  time  represents  a  transi¬ 
tion  between  two  neighboring  feedback  control  regions.  Now  by  Theorems 
3.2  -  3.4  there  are  no  returns  of  states  to  the  boundary,  there  are  no  breaks 
in  the  optimal  controls  between  junction  times  and  there  is  only  one  subregion 
per  region.  Therefore  every  trajectory  is  characterized  only  by  departures 
of  state  variables  from  the  boundary,  the  boundary  junction  times  might  be 
arbitrarily  taken  and  between  two  successive  boundary  junction  times  the  set 
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of  optimal  controls  does  not  change.  Thus,  (Operation  1  takes  into  account  a^i 
possible  trajectories  in  the  state  space,  sc  that  the  construction  c£  reedback 
control  regions  based  on  these  sequences  insures  a  complete  covering  of  the 
state  space.  Such  a  typical  sequence  looks  as  follows: 


1 l£ )  f  i ^  )  f  *  • 

,Vl},lLpl,lLp-l'’"' 

,lLf-n*l ' 

The  number  of  possible  sequences  Q(n)  is  given  by 

.Q£ri)  *  l  1  (-l)Vi(r-i)11 

r=l  i=o  v1' 

(see  [El],  p  68). 

Notice 

that  the  complexity  of  Operation  1  ±s 

of  the  order  of  nn,  so 

that  clearly  the  algorithm  is 

lmplementable  only  for 

"smali"  networks 

Example  4.1 

Consider  a  network  with  three  state  variables 

sequences  are: 

tf-l 

x^,  x2  and  Xy  The 

rf-  2 

1. 

X1 

X2 

X3 

2. 

X1 

x3 

x2 

3, 

X1 

(x2,x3) 

- 

4. 

X2 

X1 

X3 

5. 

X2 

X3 

Xi 

6. 

*2 

(x1,x3) 

- 

7. 

X3 

X1 

X2 

8. 

x3 

X2 

X1 

9. 

X3 

- 

10. 

(,x1»x2) 

X3 

- 

11. 

(Xi.Xj) 

x2 

- 

12- 

(X2 

xi 

- 

13 

^xl'x2>x3) 

- 

□  Example  4  • 
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Operation  2 

Derive  the  networks  corresponding  to  every  possible  combination  of 
state  variables  leaving  the  boundary.  Apply  the  algorithm  of  Maximal  Flow  to 
the  networks  and  find  the  corresponding  flows. 

□  Operation  2. 

From  Operation  2  we  obtain  the  first  optimal  solutions  to  the  linear 
programs  (as  explained  in  Parts  B  and  C  of  Section  2).  Also,  the  operation 
provides  the  essential  tools  for  obtaining  the  simplifications  of  the  algor¬ 
ithm  as  will  be  shown  in  Operation  3. 

The  number  of  networks  for  which  the  Maximal  Flow  Algorithm  is  to  be 
applied  is  2n-l.  The  conplexity  of  this  operation  is  therefore  exponential 
in  the  numbeT  of  nodes. 

In  Appendix  B  we  provide  a  computer  program  in  Fortran  for  finding 
the  maximal  flow  in  these  networks,  based  on  the  algorithm  of  Edmons  and 
Karp  (see  [HI]). 

Example  4.2 

All  the  possible  combinations  of  states  leaving  the  boundary  in 
example  4.1  are: 

(x^ } , {X2 } > (x^) » (x^ ,X£ } * {x^ jX^} > {X2  ^3 ^  *  ^^1  * 

□  Example  4 . 2 

Operation  3 

For  every  sequence  of  state  variables  leaving  the  boundary  listed  in 
Operation  1,  and  considering  the  results  obtained  in  Operation  2,  execute  the 
following  steps: 
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(a)  Consider  every  set  of  the  sequence  and  check  if  there  is  a  sub* 

set  of  states  B  such  that  Be  L  ,  and  such  that  there  is  a  set  B' 

P 

where  I  U  B  c  B' ,  and  such  that  the  same  value  of  maximal  flow 
P 

corresponds  to  both  sets ,  I  U  B  and  B ' . 

(b)  If  (a)  is  satisfied,  check  if  in  a  set  Lq  of  the  same  sequence, 
where  q  <  p  there  is  at  least  one  state  variable  such  that 
x.  eBV(IpJ  B). 

(c)  If  (b)  is  satisfied,  erase  from  the  list  of  Operation  1  the  sequence 


being  currently  checked. 


□  Operation  3. 


In  order  to  justify  Operation  3,  consider  first  the  following  notation: 


C£-l},{If-2>*  *  "  »{Ip}»  { 


(4,1) 


The  sequence  (4.1)  represents  all  trajectories  in  the  state  space  character¬ 
ized  by  the  sets  of  states  travelling  on  interior  arcs  in  every  time  interval. 
In  (4.1),  the  notation  (  ^  ^  (B) )  means  that  the  set  of  states  B  is  such 

that  B  C  L  ,  that  is ,  the  states  in  B  leave  the  boundary  at  t  .  Notice 
P  p 

that  Bel  , . 


Now  consider  the  following  sequence  satisfying  (a)  and  (b) : 

(If.!) . (Ip>,  ap.1(B)},...,{Iq},{(Iq.1(xj)},Uq.2>,... 

where  x^  e  B'/(IpU  B)  . 

Now  consider  the  following  sequence:* 


(4.2) 


(If.!)*-., (Ip>.  (Ip.1(B)},...,{IqUx.},{Iq.1},{Iq_2},...  .  (4.3) 


The  only  difference  between  sequences  (4.2)  and  (4.3)  is  that  in  (4.2), 

x.  leave  the  boundary  at  t  while  in  (4.3)  x.  leave  the  boundary  at  t  . . 

3  q  1  q*l 

Notice  that  the  existence  of  a  sequence  as  (4.3)  is  insured  since  Operation  1 
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requires  the  listing  of  every  possible  sequence . 

We  claim  that  all  feedback  control  regions  constructed  from  sequence 
(4.2)  are  also  constructed  from  sequence  (4.3)  so  that  sequence  (4.2)  is  redun¬ 
dant.  The  reasons  for  that  appear  in  the  following  arguments: 

By  the  considerations  of  Part  D  in  Section  2,  we  know  that  in  order  to 

find  all  optimal  solutions  (to  construct  the  corresponding  feedback  control 

region)  corresponding  to  the  time  interval  (t^ ,tg) ,  we  find  all  optimal 

solutions  corresponding  to  the  states  in  I  leaving  the  boundary,  and  among 

g 

these  solutions  we  choose  . those  that  maximize  the  expression 


JiWi  • 

i  g 

Now  it  is  not  difficult  to  see,  that  by  Theorem  3.2  the  costate  values  in  both 

sequences  (4.1)  and  (4.2)  are  the  same  at  all  times.  Moreover,  note  that  if 

all  optimal  solutions  corresponding  to  the  states  in  I^U  B  leaving  the  boundary 

are  also  optimal  solutions  to  the  problem  corresponding  to  the  states  in  B' 

leaving  the  boundary,  then  all  optimal  solutions  corresponding  to  the  states  in 

IpU  BUD  leaving  the  boundary  (for  some  set  of  states  D)  are  also  optimal 

solutions  of  the  problem  corresponding  to  the  states  in  I  V  BVDVE  leaving 

the  boundary,  where  ECB'/I  U  B.  The  above  is  easily  seen  by  noting  that  all 

P 

minimal  cuts  corresponding  to  I^U  B  are  also  minimal  cuts  of  the  network  cor¬ 
responding  to  IpUBUE,  therefore  all  minimal  cuts  corresponding  to  IpUBUD 
are  also  minimal  cuts  of  the  network  corresponding  to  IpU  BU DUE. 


Now  returning  to  sequences  (4.2)  and  (4.3),  we  have  that 

I2  -  I1  U  x.  , 

q  q  J 


(4.4) 


1  2 

(where  I*  and  1^  denote  the  sets  of  states  travelling  on  interior  arcs  in 
the  time  interval  (t^,t^)  in  sequences  (4.2)  and  (4.3)  respectively). 


Moreover  we  can  write 


ij-  aJ/ipUSK-y-'B  ,  (4.s) 

so  that  from  (4.4)  and  (4.5)  we  have 

ll  *  (Iq/V  B)U  IpUBUxj  *  (4-6) 

But  from  the  above  reasoning  all  optimal  solutions  corresponding  to  (4.5)  are 
also  optimal  solutions  of  the  problem  corresponding  to  (4.6),  moreover: 

ig-  ig  .  fell  • 

Therefore  the  sequence  (4.2)  is  redundant. 

Example  4.3 

If  in  example  4.2  we  obtain  the  same  maximal  flow  value  in  the  cases 
corresponding  to  the  sets  of  states  x^  and  {x^x^x^}  leaving  the  boundary, 
then  we  erase  from  the  list  of  example  4.1  the  sequences  numbered  (1),  (2),  (3), 
(4),  (7),  (10)  and  (11). 

□  Exanple  4.3. 

Operation  4 

For  every  state  variable  combination  leaving  the  boundary  at 
for  which  there  is  not  a  combination  L.£  such  that  Lf  C  ,  and  the  maximal 
flow  values  corresponding  to  both  combinations  are  equal,  find  all  the  optimal 
basic  solutions  of  the  corresponding  linear  programing  problem  starting  with 
an  optimal  (extended)  basic  solution  as  explained  in  Parts  B  and  C  of  Section  2, 
and  using  the  algorithm  provided  in  Part  B  of  this  section. 

□  Operation  4 . 

The  results  obtained  from  Operation  4  are  the  basic  data  needed  for  con¬ 
structing  feedback  control  regions  (the  construction  is  carried  out  in  Opera- 


For  every  remaining  sequence  (after  execution  of  Operation  3)  from  the 
list  of  Operation  it  carry  out  the  following  steps  for  every  boundary  junction 
time  tp,  in  a  sequential  order,  starting  at  t^: 

(a)  Set 

LCt)  *-1,  Vx.e  Ipml 
V  r  e  (—  ,tp) 

X .  (t)  »  0 ,  V  x.  e  B 

l  l  p-1 

*T£( tp.i»tp) 

with 

^i  (tf )  »  0,  V  {xji  e  W}  . 

(b)  Set  arbitrarily  t  -  *  1  , 

(c)  From  among  all  solutions  (y*)  obtained  in  Operation  4  with  *  I  ^ 
choose  those  that  maximize  the  expression: 


l  »icVi)yi 


X.  el  , 
1  p-1 


(d) 


Transform  the  solutions  obtained  in  (c)  to  the  set  of  rays  V  ^  (see 
[1],  p.  58)  and  construct  the  convex  polyhedral  cone: 


R  -  •  C  (R  v  V  J/R  , 
p-1  ov  p  p-1  p 

where  C  (•)  denotes  the  convex  hull, 
o 

(e)  Consider  all  sequences  that  are  identical  to  the  current  sequence 

until  the  current  time  interval  (t  ,,t  ).  The  costate  values  of 

p-1  p 

such  sequences  in  interval  (t  ,  ,t  )  will  be  identical  to  the  co- 

p-1  p 

state  values  corresponding  to  the  considered  sequence,  and  the  cor¬ 


responding  feedback  control  region  will  be  built  only  once.  This 
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insures  that  we  avoid  multiple  construction  of  the  same  control  region. 

□  Operation  5 . 

Step  (a)  follows  from  the  costate  trajectories  depicted  in  Theorem  3.2. 
Notice  that  in  the  sequences  left  after  Operation  3,  the  costate  starts 
to  change  only  from  the  boundary  junction  time  in  which  its  corresponding 
state  leaves  the  boundary,  so  that  the  calculation  of  the  costate  values  is 
an  easy  task.  Step  (b)  follows  from  the  fact  that  there  is  only  one  subregion 
per  region  with  respect  to  any  set  of  state  variables  leaving  the  boundary  (as 
proved  in  Theorem  3.4),  so  that  we  may  arbitrarily  choose  the  boundary  junction 
times.  Step  (c)  follows  from  Part  D  of  Section  2.  Step  (d)  follows  from  the 
Constructive  Dynamic  Programming  Algorithm  stated  in  [M2],  Finally,  step  (e) 
is  a  simple  conclusion  that  saves  superfluous  computational  work, 

B.  The  Method  for  Finding  all  Solutions 
to  the  Linear  Programming  Problems 

One  of  the  major  problems  in  the  implementation  of  the  algorithm 
described  in  Part  A  of  this  section  arises  from  the  need  of  finding  all  solu¬ 
tions  to  every  linear  program  required  in  Operation  4  of  the  algorithm.  In 
geometrical  terms,  the  problem  is  to  find  all  the  extremal  points  of  the  con¬ 
vex  polyhedron  formed  by  the  intersection  between  the  enlarged  Hamiltonian 
and  the  y- constraint  figure  lying  in  the  positive  orthant  of  the  y- space. 
Fortunately,  the  above  linear  programs  are  in  general  highly  degenerate  so 
that  as  it  was  clarified  in  Part  C  of  Section  2,  the  number  of  finear  programs 
we  actually  have  to  solve  is  greatly  reduced. 

The  problem  of  finding  all  solutions  to  linear  programs  has  been 
widely  investigated.  See  for  exanple  [Bl],  [Chi],  [Ch2],  [M4J,  (Rl]-  In  [B2], 
we  find  a  report  on  computational  experience  gained  using  the  algorithms 
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presented  in  [Chi]  and  [Ch2],  as  applied  to  optimal  routing  problems.  From 
this  report  it  is  seen  that  even  for  very  small  size  networks  the  amount  of 
computation  and  memory  required  is  excessive,  and  also  the  numerical  sensitivity 
of  the  algorithm  is  extremely  high- 

We  proceed  now  with  the  description  of  the  method  to  find  all  the  solu¬ 
tions  of  the  linear  programs  that  we  propose  and  later  we  shall  discuss  its 
advantages.  The  algorithm  is  based  on  the  method  of  pivoting  on  networks. 

For  a  review  of  this  method  the  reader  is  referred  to  Appendix  A. 

Starting  at  an  optimal  extended  basic  solution,  the  algorithm  picks  a 
nonbasic  link  and  looks  for  the  cycle  formed  by  basic  links  and  the  considered 
nonbasic  link.  There  always  exists  such  a  (unique)  cycle  since  every  extended 
basic  solution  corresponds  to  a  spanning  tree  in  the  network.  In  order  to 
find  the  cycle  the  algorithm  utilizes  a  labeling  technique  that  will  be 
described  later  on.  After  finding  the  cycle  we  perform  a  pivot  on  it,  reach¬ 
ing  in  this  way  either  a  new  extremal  solution  or,  if  the  cycle  is  degenerate, 
another  representation  of  the  original  extremal  solution.  If  the  degeneracy 
is  of  a  high  degree,  we  may  find  with  a  single  pivot  several  different  representa¬ 
tions  of  the  extremal  solution.  Every  new  solution  or  representation  obtained 
is  numbered  and  kept  in  memory  if  and  only  if  it  does  not  exist  there  already. 

The  above  process  is  performed  for  all  nonbasic  links  corresponding  to  every 
solution  and/or  representation  kept  in  memory  until  it  is  not  possible  to 
reach  a  new  solution  or  representation.  In  this  way  we  insure  that  at  the 
end  of  the  process  we  had  reached  all  optimal  solutions.  Since  we  are  not 
concerned  with  different  representations  of  extremal  solutions  when  construct¬ 
ing  feedback  control  regions  at  the  end  of  the  process,  we  can  delete  these 
representations  retaining  only  one  representation  per  extremal  solution.  In 
order  to  start  the  process,  we  reach  an  initial  extended  basic  optimal  solution 
as  depicted  in  Parts  B  and  C  of  Section  2. 


Now,  to  understand  the  need  of  searching  also  the  different  representa¬ 
tions  of  every  extremal  solution,  let  us  consider  the  following  exanple. 

Example  4.4 

Given  the  network  depicted  in  Figure  4.1,  the  problem  is  to  find  all 
optimal  solutions  corresponding  to  all  states  leaving  the  boundary. 


After  finding  the  minimal  cut  of  the  network  (recall  that  to  find  a 
first  optimal  basic  solution,  we  add  a  new  node  (the  source)  and  links  with 
no  capacity  constraints  connecting  the  source  with  every  node  of  the  network, 
and  we  apply  the  maximal  flow  algorithm) ,  we  choose  the  first  basic  optimal 
solution  shown  in  Figure  4.2-a.  From  this  solution  and  by  pivoting  on  the 
cycle  corresponding  to  the  nonbasic  link  (3,2),  we  reach  another  extremal 
solution  as  depicted  in  Figure  4.2-b.  Note  that  from  the  first  solution, 
pivoting  on  the  cycle  corresponding  to  the  nonbasic  link  (s,2)  does  not  lead 
to  a  new  solution  but  just  to  another  representation  of  the  first  optimal 
solution.  Clearly,  from  the  representation  of  the  extremal  solution  of 
Figure  4.2-b  we  cannot  reach  a  new  extremal  solution.  Only  starting  at  tie 
representation  shown  in  Figure  4.2-c  can  we  reach  the  new  extremal  solution 
depicted  in  Figure  4.2-d. 


-  Obtaining  extremal  points  in  the  netwotk  of  Exanple  4.4. 


Figure  4 


□  Exanple  4.4. 


Now  we  describe  the  labeling  technique  used  by  the  algorithm  to  find 
pivot  cycles.  First  we  assign  a  nuriber  to  every  node  and  link  in  the  network. 
For  exanple  see  Figure  4.3: 
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Second,  we  define  the  following  arrays: 

C (I)  =  Capacity  of  the  link  I. 

F(I,J)  *  Flow  an  the  link  I  at  the  J-th  solution. 

BASE(I,J)  ■  Type  of  the  link  I  (basic  or  nonbasic)  at  the  J-th  solution. 

LOUT  (I)  *  Exit  node  of  link  I. 

LIN(I)  =  Entrance  node  of  link  I. 
where : 

1  if  the  link  I  is  basic 

BASE(I,J)  =  ■ 

0  if  the  link  I  is  nonbasic  . 

Now  if  link  K  is  nonbasic,  in  order  to  find  the  pivot  cycle  we  search  the 
arrays  LOUT (I)  and  LIN (I)  and  we  find  all  the  incident  links  to  node 
LIN  (K)  (where  LIN(K)  is  the  entrance  node  of  the  link  K).  After  find¬ 

ing  all  these  links  we  consider  only  those  that  are  basic,  and  for  every  one 
of  them  we  carry  out  the  same  procedure  until  we  obtain  a  basic  link  entering 
(exiting)  into  (from)  node  LOUT(K).  To  every  link  searched  in  the  above 
procedure,  the  number  of  the  preceding  link  of  the  cycle  is  matched  so  that 
at  the  end  of  the  procedure  a  well-defined  cycle  is  obtained.  After  finding 
the  cycle,  and  by  checking  the  arrays  C(I)  and  F(I,J)  we  obtain  the  maximal 
change  of  flow  on  the  cycle  and  we  perform  the  pivot.  As  said  before,  when 
the  cycle  is  degenerate,  we  reach  all  possible  different  representations  by 
performing  "dummy"  pivots,  that  is  declaring  the  nonbasic  link  as  basic,  and 
the  basic  link  that  does  not  allow  to  change  the  flow  as  nonbasic. 

Standard  labeling  techniques  assign  numbers  only  to  the  nodes  of  the 
network  so  that  the  links  are  identified  by  a  two-dimensional  array  (I,J), 
where  I  is  the  exit  node  of  the  link  and  J  its  entrance  node.  However, 
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note  that  in  order  to  use  this  technique  we  must  define  in  our  case  two  arrays 
that  are  three-dimensional,  the  first  accounting  for  the  flows  on  the  links, 
the  second  for  the  type  of  links.  Moreover  a  two-dimensional  array  for  the 
link  capacities  is  needed,  thus  making  use  of  an  excessive  and  wasteful 
amount  of  memory,  especially  when  the  number  of  solutions  and  representations 
is  large.  Clearly,  by  defining  only  two-dimensional  arrays  (as  we  do)  we 
save  a  large  amount  of  memory.  Notice  also  that,  when  finding  all  the  solu¬ 
tions  to  the  linear  programs  by  means  of  pivoting  on  networks,  every  solution 
or  representation  is  identified  by  a  solution  vector  and  a  link- type  vector 
instead  of  a  tableau  (as  is  done  when  utilizing  simplex  methods)  saving  again 
a  great  deal  of  memory. 

The  main  problem  of  the  method  is  that  we  do  not  know  a  priori  the 
mmber  of  solutions  and  representations  of  the  linear  program,  so  that  we  do 
not  know  the  dimension  needed  for  the  arrays  that  will  contain  the  solutions 
or  representations.  Truly,  we  can  find  an  upper  bound  to  the  number  of  solu¬ 
tions  and  representations  since  every  basic  solution  corresponds  to  a  spanning 
tree  and  it  is  possible  to  calculate  the  number  of  spanning  trees  in  the  net¬ 
work,  however,  not  every  spanning  tree  corresponds  to  a  feasible  basic  solu¬ 
tion.  For  example,  in  the  network  depicted  in  Figure  4.4a  non-feasible  span¬ 
ning  tree  is  shown. 

nonbasic  link 
basic  link 


Figure  4.4  -  Example  of  a  non-feasible  spanning  tree. 
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Therefore,  in  general,  the  number  of  feasible  solutions  is  less  than 
the  number  of  spanning  trees,  and  for  practical  purposes  the  dimension  of  the 
arrays  having  a  coordinate  defining  the  number  of  solutions  and  representa¬ 
tions  is  estimated  and  fixed  accordingly. 

We  present  now  ‘several  examples  that  will  test  the  efficiency  of  the 
method.  All  examples  apply  to  the  case  when  all  the  states  leave  the  boundary 
(that  is,  the  case  in  which  we  obtain  the  maximal  number  of  solutions)  and 
they  were  solved  by  means  of  a  Fortran  computer  program  (see  Appendix  B) 
based  on  the  method  described  above.  All  programs  were  run  on  an  IBM  370/168 
computer,  and  for  every  example  the  maximal  memory  region  provided  was  512K. 

In  Figure  4.5  we  show  the  tested  networks  with  the  corresponding 

numbers  of  optimal  solutions  and  representations,  and  the  CPU  time  consumed. 

3 

Number  of  solutions:  24 

Number  of  solutions 
and  representations :  48 

CPU  time:  1.78  sec. 

(a) 

3 

Number  of  solutions:  20 

Nunber  of  solutions 
and  representations :  58 

CPU  time:  1-59  sec. 


Cb) 


Number  of  solutions:  60 

Number  of  solutions 
and  representations:  156 


CPU  time: 


5.41  sec. 


Figure  4.5  -  Example  of  finding  all  optimal  solutions. 


From  the  examples  of  Figure  4.5  we  obtain  an  insight  into  the  change  in  the 
number  of  feasible  spanning  trees  when  changing  capacities  of  some  links  (see 
4.5a  and  4.5b)  and  the  great  increase  in  the  number  of  feasible  spanning 
trees  when  adding  a  single  link  to  the  network  (see  4.5b  and  4.5c,  4.5d  and 


Note  that  since  the  capacities  are  integers  and  the  only  arithmetical 
operations  in  the  method  are  additions  and  subs tractions,  no  numerical  sensi¬ 
tivity  problems  arise. 

Now  we  anal  ze  briefly  the  complexity  of  the  method  for  finding  all 
optimal  solutions.  Consider  first  the  non- degenerate  case,  that  is,  the  case 
in  which  every  optimal  solution  satisfies: 

Yi  >  0,  i  ■  1,2,. ..,n  ,  (4.7) 

where  n  is  the  nunber  of  nodes  in  the  network.  Denote  also  by  i  the 
number  of  links  in  the  network.  Clearly,  if  (4.7)  is  satisfied  at  every 
optimal  solution,  the  links  corresponding  to  the  slack  variables  "y"  are  basic, 
and  the  remaining  links  are  non-basic.  Since  the  flow  on  a  non-basic  link 
(i,j)  must  be  only  zero  or  c^j ,  the  number  of  optimal  (basic)  solutions 
is  2  .  This  explains  the  great  increase  in  the  nunber  of  solutions  when 

adding  a  single  link  to  the  network.  In  the  degenerate  case  we  can  only 
estimate  an  upper  bound  for  the  nunber  of  solutions  and  representations.  This 
upper  bound  is  given  by  multiplying  the  nunber  of  spanning  trees  of  the  net¬ 
work  by  2^  The  conclusion  arising  from  the  above  arguments  is  that 
owing  to  its  exponential  complexity  the  method  can  actually  be  implemented 
only  for  relatively  small  networks. 

C.  Example  of  Feedback  Solution  to  the  Dynamic  Routing  Problem  of  a 

Single  Destination  Network  with  all  Unity  Weightings  in  the  Cost  Functional 

In  order  to  illustrate  the  algorithm  described  in  Part  A  of  this  section 
we  present  here  an  example  of  a  four- state  variable  network  for  which  we  execute 
some  of  the  operations  of  the  algorithm.  The  network  is  depicted  in  Figure  4.6. 


to 
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Here  we  calculate  the  maximal  flow  values  of  the  networks  corre¬ 
sponding  to  all  possible  combinations  of  states  leaving  the  boundary. 

This  is  carried  out  by  using  the  Maximal  Flow  computer  program  provided 
in  Appendix  B. 

In  the  case  where  all  states  leave  the  boundary,  there  is  no  need 
to  apply  the  computer  program  since,  as  stated  in  Part  B  of  Section  2,  the 
solution  is  trivial  (see  Figure  4.7). 


The  solution  is: 

“li  =  Ui6  =  Ci6  ’  1  *  2’3’4’5 

uik  81  0  ,  k  *  2, 3, 4,5 

The  value  of  the  maximal  flow  is  9. 
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The  links  belonging  to  a  minimal  cut  are  marked  by  asterisks  (•) 
each  of  the  above  programs. 


Now  we  suninarize  the  results  we  obtained: 
D 

. X2  y X3 , X4 
X1.X3.X3 

xl»x2,x4 

xl,x3,x4 

VX3 
Vx4 
X1 


x2,x3 


x2,x4 


x3’x4 


k(D) 

9 

9 

9 

9 

8 

8 

7 

7 

6 

6 

6 

5 

3 

3 


This  table  will  assist  us  in  the  execution  of  Operations  3,  4  and  5. 
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Operation  3 

In  this  operation  we  erase  redundant  sequences  from  the  list  of  Opera¬ 
tion  1.  The  Operation  is  made  by  searching  all  the  sequences  of  the  list  and 
by  checking  if  the  conditions  of  redundancy  are  satisfied  considering  the 
results  obtained  in  Operation  2. 

In  our  example  we  obtain  that  only  18  sequences  from  the  75  listed  in 
Operation  1  are  nonredundant.  These  sequences  are  marked  and  numbered  in  the 
list  of  Operation  1. 

Operation  4 

Here  we  find  all  basic  optimal  solutions  for  all  combinations  of 
states  leaving  the  boundary.  By  the  summarizing  table  of  the  maximal  flow 
values  we  realize  that  it  is  enough  to  find  all  the  solutions  for  only  five 
combinations  of  states.  Notice  that  in  the  cases  where  only  one  state  leaves 
the  boundary  the  solution  is  unique,  and  it  is  given  by  the  result  obtained 
in  Operation  2.  The  five  networks  for  which  we  have  to  find  all  optimal 
solutions  correspond  to  the  largest  combinations  of  states  (containing  at 
least  two  states)  having  different  maximal  flow  values. 

Using  the  results  of  Operation  2  (where  we  localize  minimal  cuts)  and 
by  the  considerations  of  Part  C  of  Section  2,  we  obtain  first  basic  optimal 
solutions  to  every  one  of  the  networks,  as  shown  in  Figure  4.8, 


A; 

Figure  4.8  -  The  first  optimal  solutions 


Applying  the  method  for  finding  all  the  solutions  described  in  Part  B 
of  this  Section  (the  computer  program  implementing  it  is  provided  in  Appendix  B), 
on  the  above  networks,  we  obtain  the  following  results: 
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eratiom  5 


In  this  Operation,  the  feedback  control  regions  are  constructed  with 
the  aid  of  the  results  of  Operation  4,  and  by  calculating  the  costate  values 
in  every  interval  of  every  (non- redundant)  sequence.  In  this  example  we  re¬ 
strict  ourselves  to  finding  the  costate  values  since  the  remaining  steps  do 
not  represent  interesting  features  of  the  algorithm. 

As  we  said  in  Part  A  of  Section  4,  when  developing  the  algorithm  it 
turns  out  that  the  calculation  of  the  costate  values  is  almost  trivial  for  the 
non- redundant  sequences.  In  the  following  table  we  rewrite  these  sequences 
where  under  every  interval  we  write  the  corresponding  costate  vector 
(a1,X2»^3»^4)-  We  have  taken  arbitrarily  t^  -  t  ■  1,  and  the  states  on 
every  interval  correspond  to  Ip,  that  is  states  travelling  on  interior  arcs. 


xl»x2»x3,x4  *l,x3: 

(3. 1.2. 4)  (2,0 

x1,x2,X3,X4  x1'x3 

(2. 1.3. 4)  (1,0 


X2,X3 


(1,2, 3, 4)  (0,1 


XpX2 


xl»x2 


VX4 
Cl, 0,0, 2) 

X3,X4 

(0, 0,1,2) 


(0,0, 1,2) 

VX3’X4 

(0,2, 1,1) 


x2’x3,x4 


(O', 1,1, 2) 


(0,0, 0,1) 


(0,0,0, 1) 


(0,0, 0,1) 


(0,1, 0,0) 


(0,0, 0,1) 


(0,0, 0,1) 


Section  5 


CONCLUSIONS 


We  presented  a  new  approach  for  the  construction  of  a  feedback  solu¬ 
tion  to  the  minimal  delay  dynamic  message  routing  problem  for  single  destina¬ 
tion  networks  with  all  unity  weightings  in  the  cost  functional.  The  approach 
seems  to  be  the  most  appropriate  one  for  tackling  the  problem,  since  it  fully 
exploits  the  special  structure  of  the  constraints  matrix  and  also  provides  a 
physical  meaning  to  the  problem  by  working  in  a  framework  of  networks  rather 
than  in  an  abstract  mathematical  one. 

Several  improvements  have  been  achieved  compared  to  previous  results. 
The  first  is  the  great  reduction  in  the  nuriber  of  linear  programs  to  be  solved 
by  taking  advantage  of  the  high  degree  of  degeneracy  that  in  general  character¬ 
izes  this  kind  of  problems.  The  second  is  the  derivation  of  a  suitable  method 
for  finding  all  solutions  to  the  linear  programs,  that  does  not  require  the 
application  of  simplex  techniques  to  achieve  a  first  optimal  solution  but  only 
the  application  of  a  simple  algorithm  of  maximal  flow.  Moreover,  the  method 
reaches  all  the  remaining  optimal  solutions  by  pivoting  operations  on  a  net¬ 
work,  saving  in  this  way  a  great  amount  of  computer  memory.  The  third  improve¬ 
ment  achieved  is  that  the  new  approach  provides  the  tools  for  analyzing  the 
conplexity  of  the  problem,  giving  us  an  idea  of  the  number  of  conputational 
steps  required  for  obtaining  the  feedback  solution  to  a  given  network. 

Further  research  is  required  for  Operation  1  of  the  algorithm.  Since 
not  all  sequences  listed  in  the  above  Operation  are  used  for  constructing 
feedback  control  regions ,  the  question  is  whether  a  method  can  be  provided  to 
avoid  the  listing  of  those  sequences  not  contributing  to  the  construction  of 
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the  feedback  solution.  It  also  seems  that  the  algorithm  may  be  applicable  in 
the  case  when  constant  inputs  are  present.  However,  further  research  is  needed 
to  investigate  the  influence  of  these  inputs  on  the  different  Operations  of  the 
algorithm. 

Finally,  it  also  remains  to  extend  the  application  of  the  approach  to 
the  more  general  case  of  multidestination  networks. 
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Appendix  A 

Basic  Concepts  of 

Graph  Theory,  Maximal  Flow  and  Linear  Programming  in  Networks 

1.  Basic  Definitions 

Directed  Graph:  A  set  of  nodes  N  *  {l,2,...,n>  and  a  set  of  directed 
links  L  «  { (i,j),(k, (v,w)}  connecting  pairs  of  nodes  of  N. 

Path:  A  series  of  different  nodes  in  the  graph  where  between  two 
successive  nodes  of  the  series  there  is  a  link  connecting  them. 

Cycle:  A  path  to  which  a  link  is  added  connecting  the  first  and  the 
last  node  on  it. 

Connected  Graph:  A  graph  in  which  there  is  at  least  one  path  between 
any  pair  of  nodes. 

Tree:  A  connected  graph  not  containing  cycles. 

Spanning  Tree:  A  tree  containing  all  the  nodes  of  the  graph. 

Network:  A  connected  gTaph  in  which  to  every  link  (i,j)  corresponds 
a  positive  integer  c^  called  the  capacity  of  the  link. 

2.  Maximal  Flow  on  Networks  (see  [H1],[B3]). 

We  define  two  special  nodes  of  the  network;  the  first  is  the  source 
(s)  and  the  second  the  destination  (d) . 

The  concept  flow  on  the  network  is  defined  as  follows :  A  set  of  non- 
negative  integers  u^  is  called  flew  in  the  network  if  the  following  con¬ 
straints  are  satisfied: 
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-f  if  j  *  s 

0  if  j  +  s,d 

i  if  j  -  d 


(A.l) 


0  i  Uij  ^  Cij  -  V(i,j)eL  .  (A.2) 

The  Maximal  Flow  Problem  is  defined  as: 

Maximize  f  (A.  3) 

such  that  (A.l)  and  (A.2)  are  satisfied. 

A  Cut  separating  the  nodes  s  and  d  is  denoted  by  (X,X)  and 
defined  as  follows: 


(X,X)  -  {(i,j)eL/  ieX,jeX}  , 

where  X  is  a  set  of  nodes  in  the  network  containing  the  node  s  but  not 
the  node  d,  and  X  •  N/X. 

The  Cut  Value  is 


I  c.. 

jeX  ^  • 

je* 

A  Minimal  Cut  of  the  network  is  defined  as  a  cut  having  minimal  value. 

Theorem  (The  Maximal  Flow-Minimal  Cut,  (see  [FI])) 

The  maximal  flow  in  a  network,  between  the  source  and  the  destination 
equals  the  value  of  a  minimal  cut  separating  the  source  and  the  destination. 

From  the  above  Theorem  follows  that  if  (X^,X^)  j.s  a  minimal  cut  of 
the  network,  then  at  maximal  flow  we  have 


Uik'cut 


3.  Minimal  Cost  Network  Flows  (see  [B3],[H1]).  • 

Consider  a  network  having  n  nodes  and  m  links.  To  every  node  i 

in  the  network  corresponds  an  integer  that  accounts  for  either  the 

"supply"  of  flow  to  the  node  (if  >  0),  or  the  "demand"  of  flow  from 

the  node  (if  b^  <  0).  Assume  also  that  the  total  supply  of  flow  to  the 

n 

network  equals  the  total  demand  from  it,  that  is  £  b.  ■  0.  To  every  link 

i=l  1 

(i,k)  in  the  network  corresponds  a  number  y—  0  that  expresses  the  cost 
per  unity  of  flow  on  the  link. 

The  Minimal  Cost  Network  Flow  Problem  is  the  determination  of  the 
flows  on  the  links  such  that  the  supply  of  flow  to  the  network  satisfies  the 
demand  of  flow  from  it  at  minimal  cost.  Mathematical  the  problem  is  stated 
as  follows  (the  summations  are  taken  for  existing  links) : 

n  n 

Minimize  l  l  y..u  ,  (A.  4) 

i=l  j=l  V 

such  that 

n  n 

l  l  “w  a  »  i  *  1,2, . . .  ,n  ,  (A. 4a) 

j»l  k*l  1 

°luijlcij»  V(i»j)eL  .  (A4.b) 

The  constraints  (A.  4a),  written  in  vector  form,  are  Au  *  b  . 

Now  we  investigate  the  structure  of  the  constraint  matrix  A:  Every 
row  of  A  corresponds  to  a  node  and  every  column  to  a  link,  thus  the  dimen¬ 
sions  of  the  matrix  is  (n  *  m) .  Every  column  of  A  has  only  two  non- zero 
elements,  +1  in  the  i-th  row  and  -1  in  the  j-th  row,  where  i  and  j  are 
the  exit  and  entrance  nodes  of  the  link,  respectively.  The  rank  of  the  A- 
matrix  is  (n-1).  To  see  this,  note  that  every  row  in  A  accounts  for  the 
flow  conservation  on  the  corresponding  node,  so  that  the  row  corresponding  to 
node  r  is  given  by  the  sum  of  all  the  remaining  rows  of  A  with  minus  sign. 


These  remaining  rows  are  clearly  independent  since  the  deletion  of  the  r-th 
row  from  A  leaves  at  least  one  column  having  only  one  non- zero  element. 

Now  we  shall  deal  with  the  linear  programming  problem  (A.  4),  when  the 
links  have  no  capacity  constraints,  that  is  we  require  only  >_  0.  In  a 
linear  programming  problem,  the  variables  corresponding  to  some  of  the  columns 
of  the  constraint  matrix  are  called  nonbasic  variables,  or  independent  vari¬ 
ables,  because  their  values  are  determined  arbitrarily.  The  values  of  the 
remaining  variables  (the  basic  or  dependent  variables)  are  determined  by 
means  of  the  non-basic  variables  and  the  constraint  matrix.  The  columns  cor¬ 
responding  to  the  basic  variables  form  a  square  matrix  B,  called  the  basic 
matrix,  and  if  these  columns  are  chosen  such  that  they  are  linearly  inde¬ 
pendent,  and  if  the  variables  not  associated  with  columns  of  B  are  set  equal 
to  zero,  then  the  solution  of  the  system  Bu  -  b  is  unique  (since  B  is  non¬ 
singular)  ,  and  the  above  solution  is  an  extreme  point  of  the  constraint  figure. 
In  the  case  of  problem  (A.  4),  it  is  easy  to  see  that  from  all  links  incident 
to  a  node,  the  flow  on  one  of  them  cannot  be  determined  arbitrarily  but  by  the 
flow  conservation  equation  of  the  node.  If  we  consider  all  links  on  which  the 
flows  are  determined  by  the  remaining  links  of  the  network,  these  links  form 
a  spanning  tree.  The  reason  for  that  is  that  as  we  said  before,  the  flow  on 
one  of  the  links  incident  to  a  node  is  determined  uniquely  by  the  flows  on 
the  remaining  links  incident  to  the  node  satisfying  in  this  way  the  flow  con¬ 
servation  equation  of  the  node.  Now  the  constraint  matrix  A  has  (n-1) 
linearly  independent  rows,  that  is,  there  are  (n-1)  independent  flow  con¬ 
servation  equations,  so  that  there  are  (n-1)  links  in  the  network  on  which 
the  flow  is  determined  by  the  flows  on  the  remaining  links.  These  (n-1) 
links  cannot  form  any  cycles  because  one  can  add  flows  in  the  cycle  without 
violating  the  flow  conservation  equations  of  the  nodes  in  the  cycle.  This 
would  contradict  the  fact  that  the  values  of  these  (n-1)  links  are  uniquely 
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determined.  Moreover,  to  every  node  corresponds  at  least  one  such  link,  there¬ 
fore  the  links  with  flows  representing  basic  variables  form  a  spanning  tree. 


Next  we  describe  how  to  represent  a  nonbasic  vector  (that  is,  a  column 
of  A  corresponding  to  a  nonbasic  variable)  in  terms  of  basic  vectors.  In  order 
to  do  this,  note  that  every  column  of  A  (denote  it  by  a^ )  is  of  the  form 

a. .  *  e .  -  e .  , 

-ij  -l  -j 

where  e^  and  e^  are  unity  vectors  in  En,  the  former  having  a  +1  on  the 
i-th  row  and  the  latter  having  a  +1  on  the  j-th  row.  Now,  by  choosing  (n-1) 
linearly  independent  columns  of  A  we  form  matrix  B.  Matrix  B  represents  a 
spanning  tree  in  the  network.  Now  if  we  choose  a  nonbasic  link,  say  the  link 
(v,w),  clearly  we  will  have  a  unique  path  consisting  of  basic  links  only, 
connecting  the  nodes  v  and  w,  since  the  basic  links  of  the  network  form 
a  spanning  tree.  The  above  path  together  with  the  nonbasic  link  (u,w)  form 
a  cycle  (see  Figure  Al).  By  determining  the  cycle  orientation  as  the  direc¬ 
tion  of  the  nonbasic  link  we  will  have 


Figure  Al  -  Cycle  formed  by  adding  a  nonbasic  link  to  the  spanning  tree. 
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The  method  for  representing  a  nonbasic  vector  in  terms  of  basic  vectors 
is  then  as  follows:  First,  we  determine  the  cycle  formed  by  adding  the  nonbasic 
link  to  the  spanning  tree,  Second,  we  assign  to  the  cycle  an  orientation, 
according  to  the  direction  of  the  nonbasic  link.  Third,  to  all  basic  links  in 
the  cycle  having  the  same  direction  of  the  cycle  we  assign  a  -1  coefficient  in 
the  representation  and  to  all  basic  links  in  the  cycle  with  direction  opposite 
to  that  of  the  cycle  we  assign  a  +1  coefficient  in  the  representation.  To  all 
remaining  links  in  the  spanning  tree  we  assign  a  zero  coefficient.  The  above 
coefficients  correspond  to  those  appearing  in  the  columns  of  the  Simplex 
tableau.  As  said  before,  these  coefficients  are  +1,  -1  or  zero.  This  property 
of  the  A-matrix  is  called  the  unimodularity  property  and  it  insures  that  every 
basic  solution  is  formed  by  integers  only,  provided  that  b  is  an  all- integers 
vector. 

From  the  above  considerations  it  turns  out  that  in  order  to  perform  a 
pivot  in  problem  (A.  4)  t  it  is  enough  to  foim  a  cycle  and  to  increase  the  flow 
on  every  link  of  the  cycle  in  the  direction  of  the  nonbasic  link,  until  the 
flow  on  one  of  the  basic  links  becomes  zero.  This  link  leaves  the  basis  and 
the  previous  nonbasic  link  enters  the  basis.  When  there  is  a  basic  link  in 
the  cycle,  with  direction  opposite  to  it  and  with  zero  flow,  we  say  that  the 
cycle  is  degenerate.  In  this  case  we  cannot  perform  a  pivot  changing  the 
flow  on  the  cycle,  but  we  can  obtain  another  representation  of  the  current 
basic  solution  by  removing  the  link  that  does  not  permit  flow  change  from  the 
basis  and  putting  the  nonbasic  link  into  the  basis.  If  there  are  two  or  more 
links  yielding  degeneracy,  we  can  reach  all  possible  representations  by  apply¬ 
ing  the  above  method. 

Next  we  return  to  problem  (A.  4),  but  now  consider  the  case  where  the 
links  do  have  capacity  constraints.  It  turns  out  that  the  problem  can  be 
solved  in  a  manner  very  similar  to  the  one  used  in  the  case  without  capacity 
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constraints  by  making  use  of  the  concept  "extended  basic  solution"  (see  [Ll], 
p.  48).  The  idea  is  to  treat  the  capacity  constraints  (or  upper  bound  con¬ 
straints)  in  an  inplicit  way  (similar  to  non-negative  constraints  on  variables). 
This  method  avoids  the  great  increase  of  dimensionality  in  the  problem,  result¬ 
ing  from  the  addition  of  slack  variables  to  the  upper  bound  constraints .  For 
example,  if  in  problem  (A. 4)  the  A-matrix  has  dimension  (n  *  m) ,  then  adding 
a  slack  variable  to  each  of  the  constraints  (A.  4b)  yields  to  a  constraint 
matrix  of  dimension  (n  +  m)  x  2m.  The  concept  "extended  basic  solution"  cor¬ 
responding  to  problem  (A. 4)  is  defined  as  a  feasible  solution  in  which  n 
variables  corresponding  to  linearly  independent  columns  of  A  are  basic,  and 
the  remaining  (m-n)  variables  are  nonbasic,  each  having  either  value  zero 
or  being  equal  to  its  upper  bound  (i.e.  its  capacity).  In  the  network,  an 
"extended  basic  solution"  is  characterized  by  a  spanning  tree,  that  is,  the 
graph  formed  by  the  links  corresponding  to  basic  variables  is  a  spanning  tree. 
Each  link  corresponding  to  a  nonbasic  variable  carries  either  no  flow,  or  flow 
equal  to  the  opacity  of  the  link. 

Starting  from  an  initial  "extended  basic  solution",  in  order  to  make 
a  pivot  in  the  network,  we  choose  one  of  the  links  that  is  nonbasic  and  look 
for  a  cycle  formed  by  it  with  basic  type  links.  If  the  flow  on  the  nonbasic 
link  is  zero,  then  we  check  if  we  can  increase  the  flow  on  every  link  of  the 
cycle  in  the  direction  of  the  nonbasic  link.  If  the  flow  on  the  nonbasic 
link  is  equal  to  its  capacity,  then  we  try  to  increase  the  flow  on  every  link 
of  the  cycle  in  the  direction  opposite  to  that  of  the  nonbasic  link.  In  both 
cases  we  increase  the  flow  until  either: 

(i)  The  flow  on  a  basic  link  of  the  cycle  reaches  the  value  of  zero  or  its 
capacity. 

(ii)  The  flow  cm  the  nonbasic  link  reaches  the  value  of  its  capacity 
(or  zero). 
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If  (i)  occurs  first,  then  the  basic  link  is  declared  nonbasic,  and 
the  previous  nonbasic  link  is  declared  basic.  If  (ii)  occurs  first,  then  the 
spanning  tree  does  not  change  and  the  nonbasic  link  remains  nonbasic  with 
another  flow. 

A  degenerate  situation  corresponds  to  the  case  where  there  is  at  least 
one  basic  link  in  the  cycle  with  direction  opposite  to  that  of  the  cycle  and 
with  zero  flow,  or  a  basic  link  in  the  direction  of  the  cycle,  with  flow  equal 
to  its  capacity . 
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Appendix  B 
Computer  Programs 

1.  Maximal  Flow 

The  Maximal  Flow  Algorithm  of  Editions  and  Karp  is  implemented  by  a 
Fortran  Subroutine  called  MAXFL.  The  algorithm  finds  the  shortest  path  between 
source  and  destination  on  which  an  increase  of  flow  is  feasible,  determines  the 
maximal  amount  of  flow  that  can  be  pushed  on  the  path  and  then  performs  the  in¬ 
crease  of  flow.  The  algorithm  stops  when  no  such  path  can  be  found,  and  local¬ 
izes  the  minimal  cut  that  is  nearest  to  the  source. 

Input  Data 

N  -  Number  of  nodes  in  the  network. 

C  -  Two-dimensional  array  for  the  link  capacities .  Capacity  zero 

corresponds  to  non-existing  links.  To  links  not  having  capacity 
constraints  we  assign  a  very  large  capacity  (at  least  as  the  sum 
of  the  capacities  of  the  exit  links  from  the  corresponding  node). 
The  dimension  of  C  is  (N,N). 

Program_Output 

The  program  writes  the  flow  on  every  link  corresponding  to  maximal  flow 
and  the  link  capacity.  In  addition  it  marks  with  an  asterisk  the  links  directed 
to  the  destination  that  belongs  to  the  minimal  cut,  and  writes  the  maximal  flow 
value. 

Dimens  ioned_  Local  _yariab  les 

F  *  Two-dimensional  array  for  the  link  flows.  The  dimension  is  (N,N) . 
XT,  LAB,  NODE  and  BNODE  are  one-dimensional  arrays  of  dimension  N. 
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c 

c 


1 


101 

102 

2 


103 

4 

105 

5 

10ft 

107 

3 

1  1 
12 

13 

14 
10 

1000 


MAIN  PPUGKAM 

INTEC.FR  Cift.o).F(ft,e>>  ,XTC6) 

COMMON  C.F.XT.LAB(A) .L.K.N.MFL 
DO  1000  LL*1.14 

READ <  5  ,  l  0  )  N  » ( ( C (  !  <J )  ,J=l.N).l=l  ,N) 

DO  1  I*1.N 
DO  1  J=1.N 

F( I • J)=0 
DO  2  J  — 1  •  N 
IP(C(1|J))101 .2.101 
IF(CC  J,N)>102  ,2.10? 

F<l.J)=Cf J.N) 

F<  J.N)  =C(  J.N) 

CONTINUF 
CALL  MAXFL 
WPITE16.  14  ) 

DO  3  1=1  ,N 

OO  3  js 1 ,H 

IF<C< I  103 

DO  4  I  1  =  1 .K 

1FCI-LA8( 11 ) 14.10s, 4 

CONTINUE 

GO  TO  10© 

DO  5  Jl=l.L 

IF! J-XT! J1 ) 15.107.5 

CONT 1 NUF 

WRITE!©. 11)  I .J.F{ I . J).C(  I.J) 

GO  TO  3 

WR ITE (6,12)  I  ,J.F(I .J).C(  I. J) 

CONTINUE 

WR ITE (6.13)  MFL 

FORMAT!/.9X.*l*,ll,*.*,Il.*)*.7X.I2.7X.I2) 

FORMAT !/.9X. •(•.Il,*.*.Il.*>*.7X.l?.3X.***.3X.I2> 
FORMAT  I///. 2X . • THE  VALUE  OF  THE  MAXIMAL  FLOW  IS* .14) 
FORMAT! 1 Ml »9X , *  ARC*  »  3* , 'FLOW* ,3X , *CAPACI TV* ) 

FORMAT (3© 12/36121 
CONT  I  NUF. 

STOP 

END 
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SUBROUTINE  MA*FL 

1NTEGFR  C(t  .6  >.F  <6.6  >  .XT  <  6 1  •NODE  I  6  >  .BNODE  C  6  >. DELTA 
COMMON  C.F.XT,LA?(6),L.K.N.MFL 
107  1=1 

K—Q 
L*0 

OO  7  J=1..M 
NOOEf  J  >*G 
7  BNOOE ( J ) =0 

NODE  (  1  >=<*0 
lid  DO  1  J=l.M 
Kl=NOr>E<  I  ) 

DFLT A=  IF  l  X  <  Art  f>(  FLOAT  (  K  1  )  )  >• 

IF(C( I .J) )1 00.100. 10 1 

101  !F(C(X«J) —F (  I «J>  >100*  100.10? 

102  IF(NOOF.<  J>  >  1  ,  103,1 

103  <=*♦ 1 
LAB( K >*J 

IF (C<  1  .  J>-F  (  1  ,J  ) -OF.LT  A)  105, 104. 10* 

105  DFLT A=C ( 1 , J  >—F (  I  *  J ) 

104  NOOEt J>*0ELTA 
8  NODE ( J )= I 

1F( J-Nll .121.1 
121  Jl*J 

106  IF1J1-1 1108.107. 108 
10R  I=BNODE(Jl) 

IF ( NODE ( J1 > ) 1 lO , 1 1 0, 1 09 

109  F ( I . J1 >=F( I . Jll+OELTA 

120  J1=I 

GO  TO  106 

110  Ft  J1  .1  >*Fl  J1  ,l)-OEi_TA 
GO  TO  120 

100  IF(C( J.l  111 .1 .1  1  1 

111  IFtFt J.II11  ,1  .1 12 

112  lF(NOOE(J)>l.  113.1 

113  K=K-*1 
LAH<K)=J 

IFCFt  J.l >-DFLTA> 114, 114.115 

114  OFLTA*FtJ.l J 

115  NODE  t  J )=— DELT  A 
8N00E ( J ) *  1 

1  CONTINUE 

L*L-*1 

IF (K— L> 116.117.117 
117  I=LAB(L) 

GO  TO  llfl 

116  L*0 
MFL=0 

DO  3  1=2. N 
DO  4  J  =  1  .  K 
IF (LAB ( J 1— 1 >4,3,4 

4  CONT INUE 
L=L 

XT(L  >*I 

3  CONTINUE 

DO  5  I  *1  »K 
DO  5  J*1.L 
M=CAB(  I> 

Nl*XTt  J> 

IFtCfM.Nl >>119.5.119 
119  MFL«MFL4FCM.N1 > 

5  CONTINUE 
RETURN 
END 
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2.  Determination  of  all  Optimal  Solutions 

The  algorithm  described  in  Part  B  of  Section  4  is  implemented  by  a 
Fortran  Computer  Program  conposed  of  two  subroutines  and  a  main  program. 

Genera-l  _  Flc^_  D±  agxam 


Input _  Data 

M  -  Number  of  network  links. 

C  -  One-dimensional  array  for  the  link  opacities.  The  dimension 
is  M. 


F 


Two  dimensional  array  for  the  flows  on  the  links  at  the  current 
optimal  solution.  Its  dimension  is  (M,N),  where  N  is  the 


estimated  number  of  different  solutions  and  representations.  A 
first  optimal  solution  must  be  supplied  to  the  program  (that 
obtained  from  the  maximal  flow  algorithm)  that  is  F  (1 , 1) . 

BASE  -  Two-dimensional  array  that  determines  which  links  are  basic  and 
which  nonbasic  at  the  current  solution.  Its  dimension  is  (M,N). 

If  link  I  is  basic  then  BASE  (I,*)  *  1,  if  it  is  nonbasic 
then  BASE  (I,*)  ■  0.  The  basic  configuration  of  the  first  optimal 
solution  must  be  provided  to  the  program,  that  is  BASE  (1,1). 

LOUT  -  One -dimensional  array  for  the  exit  nodes  of  the  links. 

LIN  -  One-dimensional  array  for  the  entrance  nodes  of  the  links. 

Pro  gram_  Output 

The  program  writes  all  the  extremal  optimal  solutions  and  their  number. 

In  addition  it  writes  the  nunber  of  different  solutions  and  representations 

searched  during  the  execution. 

Dimensioned  Local  yariables 

DIRF,  BFR  and  LAB  are  one -dimensional  arrays  of  dimension  M. 
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C  DETEHMINAT  ION  OF  ALL  OPTIMAL  ''OLUTIPNS  -  NA  IN  PROGRAM 

INTEGER  C(9>«F<Q,SC  J  .“ASE  ('i  .  SO  )  ,  Q  I  SF  (  «.»  >  ,  HFP  *  9  >  ,  OF.  LT 
COMMON  K  »<  1  .K2.0EL7  .  J  1  »  I  1  . C  ,LCUT  (  9  >  .  L  I N*  9  ) .F. SA SF  ,DI  R*-  . RMR. M.NB 
RE AO (5*10)  M 
READ  (  S*II)  (C(n.!-=l.\*) 

RF AO ( 5 • l  )  )  (F(! *  1  ).l *1 .V) 

RFAO(S.U)  <a4SF<  i  .  i  )  .  i  =1  ,•«  ) 

RRADC5.11>  (LOUT< I  ) ,  !  r  i  ,y j 
READCS.ll)  (LIN* I ). I =1 ,*) 

WHITE (6. 13)  <  LOUT (  I )  *u IN*  I > . ! = l • " ) 

WRITE(6.1a>  CF( I • 1) . I  =  1  .m  ) 

K  *■  I 
K1*K 
NP*I 
1 2*1 

1  J1*0 

2  J1*JI+1 

!►*  f  Jl.GT.MlGO  TO  ? 

IF  (BASE* J1 .(O.GT.OIGQ  TO  ? 

CALL.  CYCL‘.‘ 

OO  100  1*1. <2 
00  200  J-l.M 

IF  fF  (  J  *K  l  )•  NE.F  (  J*  !  )  )D  10  IDO 

200  CONTINUE 

CO  TO  2 
100  CONTINUE 

WHITE  (6. 14)  (fU,<l|,|t|,M) 

I  2=1 2+ 1 
CO  TO  2 

3  K>Ktl 

I=(<1 .GE.< )GO  TO  1 
WRITE  (6.  I*1)  12 

VRITE(6.16>  NR 

10  FORMAT ( I  2 > 

11  FORMAT <3012 ) 

1 3  FORMAT  < 1  HI . ?X  .9  <  • {  • . I  1 , • ,  • . 1  1  .  •  )  •  , 3X  )  ) 

14  FORMAT  C/.SX  .9(12.**  )  ) 

15  FORMAT C/// *2* .•NUMBER  OF  DIFFERENT  SOLUTIONS  I*. IS) 

16  F  ORMA  T  {  /  //  •  2X  «  *  NUM.3E  V  OF  OIFF*£PFNT  SOLUTIONS  AND  Rc®RF  SENTATI ONS 

is) 

DEBUG  SUBCMK 

STOP 

ENO 


Ca»«**««****t«**«******4* *  ***** *** »*****»*♦ ***  ***** ******* »* *•*«***•« 

SUBROUTINE  »IVOT 

INTEGER  C(O),F(>),S0)  .BASE  *9  •  SO  >  •  D1  »F  ( 9  I  •  PFR  (9  )  ,  OFLT 
COMMON  <»*!,* 2, OFLT  «  J1 • II  .C.LCUT* R) .LIN*  o > . F . SA^E .0  I  RE »BFR, M, NR 
K1*K1M 
00  100  1*1. M 

100  F  *  I  *K 1  >=F(  I ,K ) 

12*1  1 

1  IFCOIRF*  I  l  l.FCUl  >00  TO  2 
F*  II  ,K  1)*F(  I  1  ,«1  )-OELT 
CO  TO  3 

2  F(I1.K1>*F(11  ,K1  >4DELT 

3  IMCH.EO.JI  )CO  TO  A 
I l*OFR( I  1  ) 

GO  TO  1 

4  IMfF*  I2,K  1)  .EO.C*  12)  K,0  TO  • 

IFCF* I2.K1 ) .EO.OJGO  TO  5 

1 2  *8F  R  (  1 2  ) 

GO  TO  4 

5  OO  200  1*1. M 

200  BASE  *  I  *K  1  >*BASEf  I  .*  ) 

BASE* J1.K1)*1 
BASE* I2.K1 ) *0 
K2«K  i-I 
OO  400  Isl.K? 

OO  300  J*1.M 

IFCF  (  J.K1  ).N£.F(  J.I  )  )GO  TO  400 

IF*  BASF  <  J  »K  1  )  *NF  .!-AS‘‘  (J.l  ))G0  TO  A0C< 

300  CONTINUE 

K1*K1-1 
RETURN 

400  CONTINUE 

NR*NR4l 
RETURN 
ENO 
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»_**********»***»* ******  <««**»** ***************  ******* ****************** 

SUBPOUT l  NE  CVCUP 

1  NTEGFfi  C<Ol.F<  Q.‘  <  )  .BASF  (»  ,50  )  .  O  IRF  <  9)  *  BFP  <  O  >  ,  DE  UT  .LAB  <«*  > 

COMMON  K  *«1  ..*1.11  ,C,LOUT<  «•  )  ,LI*<  y>,F.c'ASF  .0  1  PF  •  9FR  ,M.NP 

L=0 

Ml  *1 

M2*l 

J?*Jl 

OO  100  I*1.m 
100  L^IDsO 

OI??F<  Jl  )  =  l 

1F|F  <  J1  ,K  >»NF  .C  C  J1  )  )(,•  to 
M2*? 

O 1 PF ( J 1 >=-! 

Ml  *? 

25  J*0 

1  1  =  1 

2  I«=  (  1  *GT,y  >60  TO  ft 

GO  TO  <10,1 1 >. Ml 

11  IF<LOUT<  I  >*Nt  *LOuT  (  J2  )  >60  TO  5 

GO  TO  12 

10  I F  <  LOOT  <  I  )  »  NE  .L I  N  f  .12  )  )  GO  TO  J 

12  IFfBASL<  I.Kl.NE.DoO  TO  3 
DO  200  N*1.M 
1F<LAB<N).F0. 1 JGO  TO  j 

200  CONTINUE 

DIRF< I >»1 
J*J*1 
LAB( J  )*  I 
BFR< I >=J2 
GO  TO  (2t.22)»M2 

21  IF(LIN<  I  )  .EO.LOUTf  Jl  )  >GO  TO 
GO  TO  3 

22  1F(L1N< I ) .EQ*LIN< Jl ) )GO  TO  * 

3  1*1+1 
GO  TO  2 

A  1*1 

5  IF ( I . GT »M )GO  TO  7 
GO  TO  <13 .14) .Ml 

14  IF(LIN< I ) .NE.LOUTl J? »)GO  TO  5 
60  TO  15 

13  IF<L1N<I  ».N?.LIN(  J?)  >GO  TO  6 

15  1F<BASE<  I  .A  )  1  ICO  TO  is 

OO  300  N=l,*» 

IF(LAS<N>  .EO.  I  )GO  TO  6 
300  CONTINUE 

DIPF<1  >*-l 

J*  J+l 

LA«3( J)*I 

9FP< 1 >*JL 

GO  TO  <23.24  > ,M2 

23  IF<|.OUT<  I  >,F0.L0UT<  Jl  >  )  GC»  TO  V 
GO  TO  6 

24  1F<U0UT( I ).EO*L!N< Jl ) >G0  TO  « 

6  1*1+1 

GO  TO  5 

7  L*L+1 

J2*LAB  <L  ) 

*1*1 

IF(DIPF<  J2)  .F.C.-l  >V  1*2 
GO  TO  1 

ft  I  1*  I 

0SLT«CCI>-*<  I,K> 

GO  TO  15 
«>  11*1 

DELT*F  < I  ,<  ) 

1  A  J3*9FP  < I  ) 

IF  ID  JPF<  J3  ) •£  0. !  ) GO  TO  17 
!FOEi_T*l_T.F<  J3  ,K  )  >00  TO  15 
oeLT*F<  J3.M 
GO  TO  !« 

17  I F  <  OEUT ,L  f, (C ( J? ) <  J2  , A ) >  >  GO  TO  1 n 

OELT*C< J3>-F< J3,* ) 

IP  !F< J3.IU. Jl )G0  TO  20 

I  *J3 

GO  TO  15 
CALL  PIVOT 
PETUPN 


20 


-  Ill  - 


References 


[Bl]  Balinski,  M.L. ,  "An  Algorithm  for  Finding  All  the  Vertices  of  Convex 
Polyhedral  Sets",  J.  Soc.  Indust.  Appl.  Math.,  Vol.  9,  No.  1, 
March  1961. 

[B2]  Bloom,  J.,  "A  Program  for  Chemikova's  Algorithm",  Electronic  Systems 
Laboratory,  M. I.T.,  March  1976. 

[B3]  Bazaraa,  M.S.  and  Jarvis,  J.J.,  Linear  Programming  and  Network  Flows, 
John  Wiley  §  Sons,  1977. 

[Chi]  Chemikova,  N.V. ,  "Algorithm  for  Finding  a  General  Formula  for  the  Non¬ 
negative  Solutions  of  a  System  of  Linear  Equations",  U.S.S.R.  Com¬ 
putational  Mathematics  and  Mathematical  Physics,  4,  pp.  151-158, 
1964. 

[Ch2]  Chemikova,  N.V. ,  "Algorithm  for  Finding  a  General  Formula  for  the  Non- 
negative  Solutions  of  a  System  of  Linear  Inequalities",  U.S.S.R. 
Computational  Mathematics  and  Mathematical  Physics,  5,  pp 
pp.  228-233,  1965. 

[El]  Even,  S.,  Algorithmic  Combinatorics ,  The  MacMillan  Company,  N.Y. ,  1973. 

[FI]  Ford,  L.R. ,  Jr.,  and  Fulkerson,  D.R. ,  "Maximal  Flow  Through  a  Network", 
Canadian  J.  Math. .  8  (3),  pp.  339-404,  1950. 

[HI]  Hu,  T.C. ,  Integer  Programming  and  Network  Flows,  Addison- Wesley,  1970. 

[LI]  Luenberger,  D.G. ,  Introduction  to  Linear  and  Nonlinear  Programming , 
Addison- Wes ley,  1973. 

[Ml]  Moss,  F.H. ,  "The  Application  of  Optimal  Control  Theory  to  Dynamic  Rout¬ 
ing  in  Data  Comnunication  Networks",  Ph.D.  Dissertation,  Massa¬ 
chusetts  Inst.  Technol.,  Cambridge,  1976. 

[M2]  Moss,  F.H.  and  Segall,  A.,  "An  Optimal  Control  Approach  to  Dynamic 
Routing  in  Data  Comnunication  Networks ,  Part  I :  Principles" , 

E.E.  Pub.  No.  312,  Technion  -  Israel  Inst.  Technol.,  Sept.  1977. 


■i 


-  112  - 


[M3]  Moss,  F.H.  and  Segall,  A.,  "An  Optimal  Control  Approach  to  Dynamic  Rout¬ 
ing  in  Data  Conmuni cation  Networks,  Part  II:  Geometrical  Interpreta¬ 
tion",  E.E.  Pub.  No.  319.  Technion  -  Israel  Inst.  Technol.,  Jan.  1978. 

[M4]  Mattheiss,  T.H. ,  "An  Algorithm  for  Determining  Irrelevant  Constraints 
and  all  the  Vertices  in  Systems  of  Linear  Inequalities",  Opera¬ 
tions  Research,  Vol.  21,  No.  1,  247,  Jan. -Feb.  1973. 

[Rl]  Rubin,  D.C.,  "Vertex  Generation  and  Cardinality  Linear  Programs",  Opera¬ 
tions  Research.  23,  pp.  555-564,  1975. 

[SI]  Segall,  A.,  "The  Modelling  of  Adaptive  Routing  in  Data  Conuunication 
Networks",  IEEE  Trans .  an  Conm. ,  CCM-2S,  No.  1,  pp.  85-95, 

January  1977. 


ACKNOWLEDGEMENT 


The  authors  would  like  to  thank  Dr.  F.H.  Moss  for  stimulat¬ 


ing  discussions  regarding  the  first  stages  of  this  work. 


